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INTRODUCTION 


This  report  lists  program  write-ups  for  the  multichannel 
filter  program  set  written  under  the  multichannel  filter  pro¬ 
ject  here  at  the  SDL.  All  of  these  programs  were  verified 
using  these  write-ups  which  are  new  as  free  from  errors  as 
possible.  This  set  of  programs  should  provide  a  reasonably 
complete  capability  to  analyze  signals  and  noise,  design 
multichannel  filters  to  enhance  signals  and  suppress  noise, 
evaluate  the  performance  of  these  filters  and  prepare  the 
punched  paper  tape  which  inputs  the  multichannel  filters  into 
the  Texas  Instruments  Digital  converter. 


ALPHABETICAL  LIST  OF  PROGRAMS 


PROGRAM  NAME  PURPOSE 

1.  G621  COLYTUKY  To  compute  spectra,  correlations,  and 

ordinary  coherence  functions  using  the 
Cooley-Tukey  method. 

2.  G627  COHERNCY  To  compute  and  plot  ordinary  coherences, 

phases,  and  auto  spectra  for  a  set  of 
data  channels. 

3.  G638  CONFIL  To  convolve  MCF  impulse  responses  with 

a  band-pass  filter. 

4.  G624  FKSPTRUM  To  contour-plot  frequency-wavenumber 

spectra  and  array  response  functions. 

5.  G632  HEFALUMP  To  design  and  apply  the  measured-noise 

isotropic  processor. 

6.  G640  ISOFIL  To  design  and  apply  the  theoretical 

isotropic  processor. 

7.  G639  LSTCHNCE  To  evaluate  measured-noise  and  theoretical 

isotropic  processor  and  maximum  likelihood 
MCF  performance. 

8.  Z124  MAX LI K  To  design  and  apply  the  maximum-likelihood 

MCF. 

9.  M22 0  MCFONPT  To  write  paper  tape  images  of  MCF  impulse 

response  functions  read  from  their  respec¬ 
tive  save  tapes  onto  magnetic  tape. 

10.  M2 21  MCFONPT  To  punch  the  magnetic  paper  tape  image 

written  by  program  M2  2  0  MCFONPT  on  paper 
tape  using  the  160  computer. 

11.  G630  MULTICOH  To  compute  and  plot  multiple  coherence 

functions  using  the  Cooley-Tukey  method. 

12.  G631  PARTLCOH  To  compute  and  plot  partial  coherence 

functions  using  the  Cooley-Tukey  method. 

13.  G636  PLOTFK  To  contour-plot  frequency-wavenumber 

spectra  from  an  ensemble  averaged  spectral 
matrix  read  from  tape. 

To  design  and  apply  time  prediction  and 
prediction-error  filters. 


14.  G62  8  PREDICT 


15.  G635  SPECAVE 

16.  G625  TWX 

17.  G629  VFKSPTRUM 


To  compute  an  ensemble  averaged  spectral 
matrix  from  several  samples  of  data. 

To  design  and  apply  spatial  interpolation 
and  interpolation-error  filters. 

To  contour-plot  frequtncy-wavenumber 
spectra  for  a  linear  array. 
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SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 
Title:  COLYTUKY 

COOP  Identification:  G621  COLYTUKY 
Category:  Spectral  Analysis 
Programmer:  J.  Jih 

Date:  January,  1967 

B.  PURPOSE 


This  program  computes  auto-spectra,  cross-spectra,  auto¬ 
correlations,  cross-correlations  and  coherencies  using  the 
Cooley-Tukey  method.  Fourier  transforms  of  the  data  are 
first  computed  and  then  manipulated  to  give  the  desired  re¬ 
sults.  This  program  is  intended  as  an  alternative  to  pro¬ 
gram  BLACKY. 

C.  USAGE 


1.  Operational  procedure:  This  is  a  Fortran-63  main 
program  with  the  following  subroutines:  TAPER,  SPLOT, 
SMOOTH,  DETRND,  and  COOLER.  In  addition  the  follow¬ 
ing  utility  subroutines  are  assumed  to  be  on  the 
system  tape:  SKIPREC,  ERASE,  DISC63,  and  COOL. 

2 .  Parameters : 

a)  NJOB  -  the  total  number  of  separate  job  re¬ 
quests  to  be  processed  in  this  computer 
run.  The  program  processes  job  requests 
in  sequence,  reading  in  the  required 
data  card  each  time. 


b)  ITAPE 


c)  NCODE 


the  logical  tape  number  for  the  plot 
tape.  This  must  not  be  1,  5,  or  6. 

If  this  is  0  or  blank,  no  plots  will 
be  produced,  and  the  following  two 
parameters  will  not  be  used. 

the  spacing  between  plotted  points  in 
hundreths  of  an  inch.  NCODE  =  blank, 
0,  or  1  gives  100  points/in.  NC0DE=  2 
gives  50  points/in.  etc. 


d)  IRANGE  _  the  number  of  hundredths  of  an  inch  cor¬ 
responding  to  the  maximum  value  of  the 
plots.  IRANGE  =  blank,  or  500  gives  +_ 

1  inch  plots  etc. 


e)  ISM 

f)  I CHI 


g)  ICH2 


h)  ISPT 


i)  NPTS 


the  seismogram  number  of  the  data. 

the  channel  number  for  the  first  data 
channel,  if  one  is  all  that  is  to  be 
used.  ie:  ICH1  =  4,  indicates  the  de¬ 
sired  channel  is  the  fourth  on  the  input 
tape. 

the  channel  number  for  the  other  data 
channel  if  required.  ICH2  =  0  or  blank, 
means  no  cross-spectrum,  cross-corre¬ 
lation  or  coherency  will  be  calculated. 

the  first  point  of  the  requested  data 
sample  for  each  channel  on  the  input 
tape. 

the  number  of  data  points  in  the  sample. 
If  ISPT  and  NPTS  specify  more  points 
than  exists  on  the  input  tape,  the  pro¬ 
gram  cuts  down  NPTS .  If  NPTS  is  not  a 
power  of  2,  it  is  truncated  to  the  next 
power  of  2  less  than  NPTS.  NPTS  <_  819  2. 
NPTS  £  4096,  if  correlations  are  to  be 
calculated. 


j)  ICR0SP12-switch  for  the  form  of  cross-spectra 

and  cross-correlations.  ICROSP12=0 
or  blank,  gives  ICH1  cross  ICH2. 

ICR0SP1 2  =  1  gives  ICH2  cross  ICH1. 

k)  ISMPT  -  the  number  of  points  desired  after  smooth¬ 

ing  in  the  auto  and  cross  spectra.  This 
number  is  truncated  to  the  next  power 
of  2  plus  1  less  than  ISMPT  if  it  is  not 
so  already,  ie:  17,  33,  129  etc.  If 
ISMPT  =  NPTS+1,  0  or  blank,  no  smooth¬ 
ing  will  be  done. 

l)  LAGS  -  the  number  of  lags  desired  in  the  cor¬ 

relation  function  plots.  LAGS/ 2  points 
are  plotted  on  either  side  of  the  peak. 
Correlation  functions  are  only  plotted 
in  this  program,  they  are  not  printed. 

If  LAGS  =  0  or  blank,  no  correlation 
functions  will  be  computed. 
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m)  DEMAG 1  -  demagnif icati  n  factor  for  the  first 

channel,  (I CHI).  This  is  divided  into 
every  point  of  ICH1.  If  DEMAG 1  =  Q.O 
or  blank,  it  is  set  to  1.0. 

the  fraction  of  the  record  to  apply  a 

cosine  taper  to  the  beginning.  TPS  = 

0.0  or  blank  for  no  taper. 

the  fraction  of  the  record  to  apply  a 

cosine  taper  to  the  end.  TPE  =  0.0 

or  blank  for  no  taper. 

-  switch  for  detrending  the  data.  IDT  = 

-1  for  no  detrending.  IDT  =  0  to  remove 
the  mean.  IDT  =  1  to  remove  the  mean 
and  linear  trend. 

-  switch  for  spectra  plots.  ILOG  =  0 
or  blank  gives  linear  plots.  ILOG  = 

1  gives  log  to  the  base  10  plots. 

r)  IPLOT  -  switch  for  plots.  IPLOT  =  0  or  blank 

means  no  plot  is  desired  for  this  par¬ 
ticular  request.  IPLOT  =  1  means  plot 
are  desired  for  this  request. 

s)  ICOH  -  switch  for  coherency.  ICOH  =  0  or 

blank,  means  no  coherency  is  desired  for 
this  request.  ICOH  =  1  means  coherency 
is  to  be  computed  for  this  request. 

t)  IPH  -  switch  for  phase.  IPH  =  0  or  blank 

means  no  phase  is  desired  for  this  re¬ 
quest.  IPH  =  1  means  phase  is  to  be 
computed  for  this  request. 

u)  DEMAG 2  -  demagnification  factor  for  the  second 

channel,  (ICH2). 

3.  Space  required:  16510 

4.  Temporary  space:  Disc 

5.  Alarms :  Three  alarm  messages  are  given  for  the  fol¬ 
lowing  errors: 

a)  The  requested  seismogram  number  is  not  on  the  in¬ 
put  tape. 

b)  A  requested  channel  was  not  in  the  requested  seis¬ 
mogram. 

c)  A  requested  channel  was  the  timing  trace. 
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n)  TPS 

o)  TPE 

p)  IDT 

q)  ILOG 


6. 


Error  returns:  none 


7.  Error  stops:  There  are  no  stops.  The  program  always 
proceeds  with  the  next  job  request. 

8.  Tape  mountings:  An  input  tape  of  SUBSET  seismograms 
must  be  on  logical  unit  1.  If  plots  are  desired,  a 
plot  tape  must  be  on  the  unit  specified  by  parameter 
(b). 

9.  Formats :  The  first  card  of  the  data  deck  lists  para¬ 
meters  (a)  through  (d)  as  FORMAT  (18,  213,  15).  It 

is  followed  by  NJOB  cards,  one  for  each  request,  list¬ 
ing  parameters  (e)  through  (u)  as  FORMAT  (18,  213, 

215,  12,  215,  3F5.2,  513,  F5.2). 

10.  Selective  jumps:  none 

11.  Timing:  Fast! 

12.  Accuracy :  Single  precision. 

13.  Caution  to  user:  none 

14.  Equipment  configuration:  Standard  COOP  with  E/5  = 

50/6  =  51. 

15 .  Reference: 

D.  W.  McCowan ,  Finite  Fourier  Transform  Theory  and  Its 
Application  to  the  Computation  of  Convolution,  Cor- . 
relation  and  Spectra,  UED  Research  Department  Technical 
Memorandum  Number  8-66,  December  16,  1966. 

D.  METHOD 


See  reference. 


SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 

Title:  Coherency 

COOP  Identification:  G627  COHERNCY 
Category:  Time  Series  Analysis 

Programmer:  V.  Bruffey 

Date:  March,  1967 

B.  PURPOSE 

To  compute  all  the  ordinary  coherences,  auto-spectra, 
and  phases  of  a  set  of  input  data  channels. 

C.  USAGE 

1.  Operational  procedure: 

This  program  is  written  in  Fortran-63  and  uses  the 
following  subroutines: 


a) 

MAKESUB  - 

this  subroutine  reads  an  input  seis¬ 
mogram  and  produces  another  seismo¬ 
gram  with  the  requested  data. 

b) 

XMNG 

computes  spectral-matrix. 

c) 

SMOOTH  - 

smooths  spectra  to  requested  length. 

d) 

DETRND  - 

removes  mean  and;or  linear  trend. 

e ) 

MATRA66  - 

computes  matrix  transpose. 

f ) 

COOLER  - 

computes  Fourier  series  expansion  of 
a  real-valued  time  series. 

g) 

XPLOT 

writes  plot  tape  with  scale  factor  as 
input . 

Parameters : 

a) 

ISEIS  - 

the  seismogram  number  of  the  data. 
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b )  NCH 

c)  NPTS  - 

d)  ISPT  - 

e)  IDT 

f )  LF 


g)  IPRINT  - 

h)  ICH 


i)  SF 

j )  ISPACE  - 


the  number  of  channels  to  be  processed. 

the  number  of  data  points  in  the  sample 
NPTS  <  4096.  Also  NPTS  must  be  a  power 
of  2. 

the  first  point  of  the  requested  data 
sample  for  each  channel  on  the  input 
tape . 

switch  for  detrending  the  data.  IDT  = 
-1  for  no  detrending,  IDT  =  0  to  re¬ 
move  mean,  IDT  =  1  to  remove  the  mean 
and  linear  trend. 

number  of  points  desired  after  smooth¬ 
ing.  LF  must  be  a  power  of  2  plus  1, 
and  LF  <  NPTS  +1.  If  LF  =  NPTS  +  1, 
no  smoothing  is  performed.  LF  is  trun¬ 
cated  to  next  power  of  2  plus  1  if  it 
is  not  already  so. 

switch  for  matrix  print.  If  no  matrix 
is  desired,  IPRINT  *  0.  If  IPRINT  =  0, 
or  left  blank,  the  matrix  printout  is 
given. 

an  array  of  length  NCH  which  specifies 
the  data  channel  indicies.  For  example 
if  the  input  tape  contains  channels  Zl, 
Z2,  Z6»  Z7,  and  Z9,  respectively,  and 
you  want  to  process  Zl,  Z6,  and  Z9, 
then  ICH(l)  =  1,  ICH(2)  =  3  and  IC.H(3) 
5. 

an  array  of  length  NCH  which  specifies 
the  demagnification  factor  for  each 
channel.  These  are  divided  into  every 
point  of  the  data.  If  SF(I)  =  0.0  or 
blank,  it  is  set  to  1.0. 

the  number  of  hundredths  of  an  inch  for 
each  points  along  the  frequency  axis 
of  the  plot. 


3.  Storage  required:  22644 

4.  Temporary  storage :  Disc 

5 .  Error  printouts : 

If  one  of  the  following  errors  occur,  the  program  will 
skip  that  set  of  input  data,  and  proceed  to  the  next 
set  of  inputs . 
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a)  Seismograms  not  on  tape. 

b)  Number  of  points  not  a  power  of  two,  or  number 
of  points  is  greater  than  4096. 


c)  Data  array  larger  than  20000.  (i.e., 
20000). 

NCH*NCH*LF 

d)  Too  many  channels  requested  from  seismogram. 

e)  Too  many  points  requested  from  seismogram. 

6. 

Tape  mounting: 

unit  3 . 
scratch 

An  input  tape  (subset  tape)  must  be  mounted  on  logic. 

Logical  unit  2  is  used  along  with  56  and  57  as  a 
tape,  and  logical  unit  4  is  used  for  the  plot  tape. 

7. 

Input  format: 

Card  1:  FORMAT (815) 

Columns 

Variables 

1  - 

5  (right  adjusted  integer) 

ISEIS 

6  - 

10  (right  adjusted  integer) 

NCH 

11  - 

■  15  (right  adjusted  integer) 

NPTS 

16  - 

■  20  (right  adjusted  integer) 

ISPT 

21  - 

■  25  (right  adjusted  integer) 

IDT 

26  - 

•  30  (right  adjusted  integer) 

LF 

31  - 

■  35 

IPRINT 

36  - 

40 

I SPACE 

Card  ( s  )  2:  F0RMATQ6I5) 

Columns 

Variables 

1  - 

5  (right  adjusted  integer) 

ICH(l) 

6  - 

• 

10  (right  adjusted  integer) 

I CH ( 2 ) 

• 

• 

• 

76  - 

80  (right  adjusted  integer) 

• 

ICH(16) 
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If  the  number  of  channels  is  greater  than  16,  then  a 
second  card  must  follow  with  columns  (1-5)  containing  ICH(17), 
columns  (6-10)  containing  ICH(18),  etc. 


Card ( s )  3  : 


FORMAT  (3F10.5) 


Columns 


Variables 


1-10 


SF(1) 


11  -  20 


SF(  2 ) 


71  -  80 


SF( 8 ) 


If  number  of  channels  is  greater  than  8,  then  a  second 
(and  maybe  3rd  or  4th)  card  is  used  with  (1-10)  containing 
SF(9),  11-20  containing  SF(10)  etc. 

Note:  The  program  works  for  multiple  input  sets,  but  assumes 

all  seismograms  being  processed  on  one  run,  are  on  the 
same  subset  tape.  Also,  a  single  seismogram  may  be 
processed  over  and  over  on  the  same  computer  run. 

8.  Selective  jumps:  none 

9.  Timing :  10  channels  of  4096  points  each  and  smoothing 

to  129  points  requires  about  20  minutes. 

10.  Accuracy:  single  precision 

11.  Caution  to  usdr:  NCHANCH*LF  <  20000 

LF<NPTS<  4096,  LF  =2X  +  1,  NPTS  =  2y 
for  some  x  and  y. 


Output :  The  output  is: 

(1)  LF  square  matrices  of  dimension  NCH*NCH.  The 
upper  half  of  the  matrices  contains  coherences, 
the  lower  half  contains  the  phases  scaled  by  tt 
and  the  diagonal  contains  the  auto-spectra  scaled 
by  their  largest  value,  ro  that  they  lie  between 

0  and  1.  The  print  out  is  omitted  if  IPRINT  t  0. 
If  NCH  >  16,  a  two  page  print  is  used. 

(2)  NCH  plots,  where  each  plot  has  LF*NCH  points  and 
each  plot  can  be  broken  down  into  NCH  segments. 

For  example,  if  NCH  =  3  then  plot  #1  would  contain 
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the  auto-spectra  of  channel  1,  the  coherency 
between  channels  1  and  2,  and  the  coherency 
between  channels  1  and  3.  The  second  plot 
would  contain  the  phase  between  channels  2 
and  1,  the  auto-spectra  of  channel  2,  and  the 
coherency  between  2  and  3.  Finally,  the  third 
plot  would  contain  the  phase  between  channel 
3  and  1,  the  phase  between  channels  3  and  2,  and 
the  auto-spectra  of  channel  3.  The  resulting 
plot  matrix  is  shown  below. 


Plot  1 
Plot  2 
Plot  3 


Auto-spectra (1,1) 
Phase  (2,1) 

Phase  (3,  1) 


Coherency  (1,2) 
Auto-spectra  (2,2) 
Phase  (3,2) 


Coherency  (1,3) 
Coherency  (2,3) 
Auto-spectra  (3,3) 


?,hus’  uppe^  1?alf  of  the  Plot  matrix  contains  the  coherences 
the  lower  half  contains  the  phases,  and  the  diagonal  contains 
the  auto-spectra. 


13 •  Equipment  configuration:  standard  COOP. 
D.  METHOD 


procedure°mPUtati0nal  meth°d  is  Siven  by  the  following  step 

1)  Select  subsetted  data. 

2)  Compute  the  spectral  matrix  of  the  selected  data.  This 
is  performed  by  subroutine  XMNG.  The  results  are  stored 
in  matrix  form  with  imaginary  components  above  the  di¬ 
agonal,  real  components  below  the  diagonal,  and  the 
auto-spectra  on  the  diagonal. 

3)  Compute  coherences  using 


(w) 


M 


S,j<w)  S  .  .  (w)  ’  1  *  ^ 


11 


13 


the 


Since  C . . (w) 
matrix  us 


=  C.^(w)  and  C.. (w)  =  1,  the  upper  half  of 
ed  for  storing  the  phase  relations  defined  by: 


/Imaginary  S^ (w) \ 

phase  ij  (w)  =  i  TAN-1  (  - — - L_  ...  )  £ 

it  \ Real  S.  .  (w)  /»  1 


i  /  j 


13 


4.3he^Ut°~S?eCtra  are  scaled  by  their  largest  value  and  left 
on  the  diagonal. 
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SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computing  Section 

A.  IDENTIFICATION 

Title:  Band-Pass  Filter  Convolution  Program 

COOP  Identification:  G638  CONFIL 
Category:  Time  Series  Analysis 

Programmer:  D.  W.  McCowan 

Date:  November,  1967 

B.  PURPOSE 

This  program  reads  band-pass  filter  coefficients  from 
cards  and  convolves  them  with  the  impulse  response  of  a 
multichannel  filter  read  from  tape.  It  is  designed  to  pro¬ 
cess  two  kinds  of  multichannel  filters,  the  measured-noise 
and  isotropic  processors,  from  their  save  tapes.  Any  band¬ 
pass  filter  may  be  used  provided  the  number  of  coefficients 
in  its  impulse  response  does  not  exceed  100. 

C.  USAGE 

1.  Operational  Procedure:  This  is  a  Fortran-63  main 
program. 

2 .  Parameters : 

a)  JWORD  -  a  word  indicating  the  type  of  save 

tape  to  process 

JW0RD=  6HIS0FIL  for  theoretical  MCF 
JW0RD= 8HHEFALUMP  for  measured-noise  MCF. 


b)  NFP  -  the  number  of  coefficients  in  the 

band-pass  filter  impulse  response. 
NFP<100 

c)  FIL  (I)  -  the  Ith  band-pass  filter  coefficient 

3.  Space  required:  32  K 

4.  Temporary  storage:  None 


J 


5.  Alarms .  Various  alarms  for  illegal  data  are  printed 
out 


6.  Error  returns :  None 

7.  Error  stops:  A  stop  occurs  after  each  alarm 

8‘  :  Tape  unit  1  is  the  input  save  tape  on 

which  is  stored  the  MCF  filter,  Tape  unit  2  is  the 
output  tape  in  the  same  format  as  the  input  save  tape 
but  with  the  MCF  now  convolved  with  the  B-P  filter. 

9*  Input  and  output  formats:  Parameter  (a)  appears  on 
the  first  card,  as  FORMAT  (A8).  Parameter  (b)  appears 
on  the  next  card  as  FORMAT  (110).  Parameters  (c) 
appear  on  the  following  cards  as  FORMAT  (5x,E15.7). 

10-  Selective  jumps:  None 

11.  Timing:  Fast 

12.  Accuracy :  Single  precision 

13.  Cautions  to  user:  None 

14 *  Equipment  configuration:  Standard  COOP  with  RELOCOM, 
and  two ( 2 ) tape  units . 

15.  References :  None 
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SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 
Title:  FKSPTRUM 

COOP  Identification:  G624  FKSPTRUM 
Category:  Time  Series  Analysis 

Programmer:  J.  Jih 
Date:  January  1967 

B.  PURPOSE 


To  compute  and  display  the  frequency-wave  number  power 
spectra  of  seismic  noise  along  with  a  response  function  for 
the  corresponding  array.  This  program  is  intended  as  a  re¬ 
placement  for  certain  parts  of  program,  PEAKAY  using  the 
Cooley-Tukey  method  to  estimate  auto  and  cross  power  spectral 
density  functions. 

C.  USAGE 


1.  Operational  procedure:  This  is  a  Fortran- 6 3  main 
program  with  the  following  subroutines:  COOLER, 
SMOOTH,  DETRND ,  LOCATION,  FKMAT,  and  CNTUR5 .  In 
addition  the  following  utility  subroutines  are 
assumed  to  be  on  the  system  tape:  SKIPREC,  ERASE, 
DISC63 ,  and  COOL. 


2.  Parameters: 


a)  NJOB 

b)  ISM 

c)  N 

d)  NF 


the  total  number  of  separate  job  re¬ 
quests  to  be  processed  in  this  com¬ 
puter  run.  The  program  processes  job 
requests  in  sequence,  reading  in  the 
required  data  cards  each  time. 

the  desired  seismogram  number  for  this 

job . 

the  desired  number  of  input  data  chan¬ 
nels  (<_  27). 

the  number  of  frequencies  for  which  f-k 
spectra  are  to  be  computed  (<_  24). 


e)  ISPT 


f )  LX 


g )  ISMPT 


h)  IRESP 


i)  NKX 
and 
NKY 


j )  ICH 
(I) 


k)  XAXi'S 
(I),S 
YAXIS 
(I) 

l)  SF(I) 


the  first  point  of  the  requested  data 
sample  for  each  channel  on  the  input 
tape. 

the  number  of  data  points  in  the  sample. 

If  ISPT  and  LX  specify  more  data  points 
than  exist  on  the  input  tape,  the  pro¬ 
gram  goes  on  to  the  next  job  request. 

LX  is  always  truncated  to  the  next  power 
of  2  less  than  LX  if  it  is  not  so  already. 
In  the  following  discussion  this  value 
will  be  referred  to  as  LXT„  LXT  <  4096. 

the  number  of  points  desired  after  smo¬ 
othing  in  the  auto  and  cross  power  spectral 
density  functions.  This  number  is  trun¬ 
cated  to  the  next  power  of  2  plus  1  less 
than  ISMPT  if  it  is  not  so  already,  i.e., 
17,  33,  129,  etc.  In  the  following  dis¬ 
cussion  this  value  will  be  referred  to 
as  ISMPTR .  If  ISMPT  =  0  or  blank,  no 
smoothing  will  be  done  on  the  spectra 
and  ISMPTR  will  be  set  to  LXT+1. 
ISMPTR<LXT+1<4097  . 

switch  for  requesting  an  array  response. 
IRESP>0  gives  a  response  plot  with  the 
FK  spectra. 

IRESP=  0  gives  no  response  plot  with  the 
FK  spectra. 

IRESP<0  gives  a  response  plot  only  and 
no  FK  spectra. 

the  number  of  KX  and  KY  values  res¬ 
pectively  for  which  the  f-k  spectra  and 
the  array  response  function  are  to  be 
computed.  (NXK  _<  63,  NKY  <_  63).  These 
numbers  must  also  be  odd. 

the  position  of  the  Ith  data  channel 
on  the  input  tape,  ie:  ICH(I)  =  4 
indicates  the  Ith  channel  is  the  fourth 
on  the  input  tape.  (I  <_  N). 

the  X  and  Y  coordiantes  the  Ith  channel. 

The  computed  spectra  and  response  func¬ 
tion  will  be  for  arguments  of  one  over 
these  units.  (I  <_  N). 

a  demagnifiration  or  scale  factor  for  the 
Ith  channel,  which  is  divided  into  the 
data.  (I  <  N  . 
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m) 

FREQ 

(I) 

n) 

VELR 

o) 

FKX 

and 

FKY 

P) 

DKX 

and 

DKY 

the  Ith  requested  frequency.  These 
numbers  must  be  in  increasing  numeric 
order. 

the  velocity  range  for  all  spectra  plots. 

the  lowest  values  of  KX  and  KY  for  the 
response  plot  if  requested. 

the  increment  of  KX  and  KY  for  the  re¬ 
sponse  plot  if  requested. 


3.  Space  required:  20807 

4.  Temporary  space:  Disc 

5.  Alarms :  Four  alarm  messages  are  given  for  the  follow- 
ing  errors : 

a)  The  requested  seismogram  number  is  not  on  the  in¬ 
put  tape. 

b)  The  requested  data  exceeds  the  data  on  the  input 
tape. 

c)  The  number  of  channels  exceeds  the  number  of  chan¬ 
nels  on  the  input  tape. 

d)  A  requested  channel  was  the  timing  trace. 

6.  Error  returns:  none 

7.  Error  stops:  There  are  no  stops.  The  program  always 
proceeds  with  the  next  job  request. 

8.  Tape  mountings:  An  input  tape  of  SUBSET  seismograms 
must  be  on  logical  unit  1. 

9.  Formats:  The  first  card  of  the  data  deck  lists  para¬ 
meter  (a)  as  FORMAT  (18).  It  is  followed  by  NJOB 
request  decks.  The  first  card  of  each  request  deck 
lists  parameters  (b)  through  (i)  as  FORMAT  (18,  214, 

12,  316,  215).  The  next  cards  list  parameters  (j), 

(k),  and  (1)  as  FORMAT  (3(13),  F8.2,  F6.2)).  The 
next  cards  list  parameter  (m)  as  FORMAT  (12F6.2). 

The  last  card  in  each  request  deck  lists  parameters 
(n)  through  (p)  as  FORMAT  (5F6.2).  There  must  be  no 
blanK  cards.  The  next  request  deck  follows,  if  1RESP<0 
parameters  (b) , (d) , (d) , (f ) , (g) , (1) , (m) ,  and  (n)  are  not 
needed.  Furthermore  the  card  listing  parameter  (m) 

is  not  read. 
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10.  Selective  jumps:  none 


11.  Timing :  A  NKX  =  NKY  =  21,  N  =  10,  LX  =  4096,  ISMPT  = 
257,  IRESP  =  1,  NF  =  6  case  takes  24  minutes  of  CDC 
1604B  time. 


12.  Accuracy:  Single  precision 


13.  Caution  to  user:  All  data  is  detrended 


Equipment  configuration:  Standard  COOP  with  E/5  = 
50/6 =  51 .  A  RELCOM  card  is  required. 


15.  References 


D.  W.  McCowan,  Finite  Fourier  Transform  Theory  and  Its 
Application  to  The  Computation  o^B^hvolutions ,  Cor- 
relations,  and  Spectra,  UED  Research  Department,  Tech¬ 
nical  Memorandum  Number  8-66,  December  16,  1966. 


D.  METHOD 


The  frequency-wave  number  power  spectra  are  computed  from 
the  spectral  matrix  elements  by  the  following  relation: 

N  N 


P(f,kx,ky)  =  S..(f)  exp  27ri(kx(x.-x.)  +  ky(y±-y.)) 


The  spectral  matrix  elements  are  computed  as  described  in 
the  reference. 
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SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 


Title:  Measured  Noise  Isotropic  Processor 

COOP  Identification:  G632  HEFALUMP 
Category:  Time  Series  Analysis 
Programmer:  D.  W.  McCowan 

Date:  September,  1967 

B.  PURPOSE 


This  program  computes  and/or  applies  a  multichannel  iso¬ 
tropic  processor  to  seismic  array  data.  An  actual  noise 
model  is  used  computed  from  the  spectra  of  a  specified  data 
sample.  Either  a  point  or  a  disc  signal  model  can  be  com¬ 
puted.  The  program  then  solves  the  multichannel  Wiener-Hopf 
equation  in  the  frequency  domain  to  get  the  optimum  filter 
which  rejects  the  noise  and  passes  the  signal.  The  filter 
is  written  on  a  save  tape  for  future  use.  An  option  exists 
to  filter  a  given  piece  of  data  and  plot  the  filtered  trace 
and  the  direct  sum.  PERSON  plots  are  also  generated  and 
plotted  which,  in  conjunction  with  the  output  plots,  allow 
the  user  to  calculate  signal  to  noise  ratios. 

C.  USAGE 

1.  Operational  procedure:  This  is  a  Fortran- 6 3  main 
program  with  the  following  subroutines:  COMFIL,  RES- 
CALAT,  COOLER,  DETRND,  MATMPY 6 3 ,  MATRA6  3 ,  SMOOTH, 

SMNG,  TRANSIT,  FOOBACK,  MAKESUB,  BES,  COOLCON,  PLOTEM, 
PERS0N66 ,  and  SPLOT.  In  addition,  the  following 
utility  subroutines  are  assumed  to  be  on  the  system 
tape:  ERASE,  COOL,  and  DISC63. 

2 .  Parameters : 


a) 

ISEIS  - 

The  seismogram  number 

of  the  noise 

sample . 

b) 

NCH 

The  number  of  channels 

in  the  filter 

NCH  <  13. 
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c)  NPTS 


d)  ISPT 

e)  IDT 


f )  L 


g)  ISW 


h)  IND 


i)  SNR 


j )  VEL 


k)  IPER 


1)  IDISCSIG  - 


The  number  of  points  from  seismogram 
ISEIS  to  use.  This  must  be  a  power  of 
two.  NPTS  <  4096. 

The  first  point  of  the  noise  sample. 

The  detrend  switch  for  the  noise  sample 
IDT  -  -1  for  no  detrending 
IDT  =  0  to  remove  the  mean  only 
IDT  =  1  to  remove  the  mean  and  linear 
trend . 

The  number  of  points  in  the  frequency 
representation  of  the  filters.  This 
must  be  a  power  of  two  plus  one;  i.e, 

17,  33,  129  etc.  L  <_  129.  The  number 
of  points  in  the  time  representation 
of  the  filters,  the  impulse  response 
will  be  2*(L-1)-1. 

A  control  switch 
=1  to  compute  the  spectral  matrix, 
compute  the  filter  and  filter  the 
data. 

=2  to  read  the  spectral  matrix  from  the 
save  tape,  compute  the  filter  and 
filter  the  data. 

=3  to  read  the  spectral  matrix  and  the 
filter  from  the  save  tape  and  filter 
the  data. 

The  data  channel  to  use  as  a  signal 
power  spectrum. 

The  desired  signal  to  noise  ratio.  We 
have  found  SNR  =  4.0  to  work  well. 

The  signal  velocity  for  a  disc  model 
if  one  is  used.  The  filters  will  pass 
"signals'1  with  any  velocity  greater 
than  VSG. 

The  percentage  of  white  noise  to  add 
to  the  signal  model.  The  white  noise 
stabilizes  the  computation  and  assures 
good  convergence  in  the  filters.  We 
have  found  that  two  or  three  percent 
is  usually  adequate. 

A  control  switch  for  the  signal  model. 

=  0  for  an  infinite  velocity  signal  model 
=  1  for  a  disc  signal  model  V  >  VSG. 


E 


m)  ICH(I) 


n)  X(I) 


An  array  giving  the  positions  of  the 
NCH  data  channels  used  as  a  noise 
sample  on  the  input  tape.  ICH(I) 

=  4  indicates  the  Ith  channel  is 
the  fourth  on  the  input  tape.  The 
entries  in  ICH  must  be  in  increas¬ 
ing  numeric  order. 

The  X  coordinate  of  the  Ith  data 
channel. 


o)  Y(I ) 


The  Y  coordinate  of  the  Ith  data 
channel . 


p)  SF(I) 

q)  ISEIS1 


A  scale  or  demagnification  factor 
for  the  Ith  data  channel.  If  SF(I) 

=  0.0  or  blank  it  is  set  to  1.0. 

The  seismogram  number  of  the  data  to 
be  filtered. 


r)  NPTS1 

s)  ISPT1 

t)  I DTI 


The. number  of  points  to  be  filtered. 
NPTS1  <  3500. 

The  first  point  to  be  filtered. 

The  detrend  switch  of  the  data  to  be 
filtered. 


u)  ICHKI) 


The  position  of  the  Ith  data  channel 
to  be  filtered. 


v)  SFKI) 


The  demagnification  factor  for  the 
Ith  data  channel  to  be  filtered. 


3.  Space  required:  32K  locations 

4.  Temporary  storage:  Disc 

5.  Alarms :  Various  alarms  for  illegal  data  are  printed 

out . 


6.  Error  returns :  none 

7.  Error  stops :  A  stop  occurs  after  each  alarm. 

8.  Tape  mountings :  Tape  unit  1  is  the  input  subset  tape 
with  both  the  'noise  sample  and  the  data  to  be  filtered. 
Tape  2  is  the  plot  tape.  Tape  unit  3  is  the  save  tape 
which  should  be  a  scratch  tape  if  ISW  =  1  or  2  and 

an  input  tape  when  ISW  =  3.  Tapes  4  and  5  are  scratch 
tapes . 
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9*  Input  and  output  formats:  Input  cards: 

Parameters  (a)  through  (1)  appear  on  the  first  card 

7I5’  E15-3-  2I3>-  Parameters 

Km)  through  (p)  appear  on  the  next  cards  as  FORMAT 
(3(13,  2F8.2,  F6.2)).  Parameters  (q)  through  (t) 
appear  on  the  next  card  as  FORMAT  (4110).  Para- 

Tsa?8  16X  ^6  2)  j  appSar  °n  the  last  cards  as  FORMAT 
Save  tape: 


The  save  tape  consists  of  one  file  of  binary  data 
written  in  three  logical  records.  The  first  record 
is  a  127  word  label  written  as  follows: 


LAB ( 1 )  = 
LAB(2)  = 
LAB ( 3 )  = 
LAB (100 ) 
LAB (101) 
LAB(102) 


ISEIS 

LAB(103)  =  IND 

LAB(104)  =  SNR 

NCH 

LAB(105)  =  VEL 

NPTS 

LAB (10 6)  =  IPER 

=  IDT 

LAB(107 )  =  IDISCSI6 

“  L 

LAB (58 — 58+NCH-l )  = 

=  ISW 

LAB ( 79--79+NCH-1)  = 

The  second  record  is  a  21801  word  array  with  the 
spectral  matrix  stored  close  packed  as  follows: 


S^j (w) 


S(i,j,w),  ( (S(i, j ,w),  i=l ,NCH) ,  j=l,NCH) ,w=l,L) 

The  third  record  is  a  5397  word  array  with  the  filter 
stored  close  packed  as  follows: 


Fi(t)  =  F(i, t) ,  ( (F(i , t) ,  i=l,NCH),  t=l, 2*(L-1)-1) 

10.  Selective  jumps:  none 

11.  Timing:  A  13  channel  129  point  filter  computed  and 
applied  to  3500  points  of  data  required  about  45 
minutes  of  CDC  1604B  time. 


12,  Accuracy :  Single  precision. 

13 *  Cautions  to  user:  It  is  not  possible  to  generate  a 
plot  of  the  correlation  matrix  with  this  program. 

14 *  Equipment  configuration:  Standard  COOP  with  RELOCOM 

15 .  References : .  Seismometer  Array  and  Data  Processing 
System,  Project  Vela-Unif orm-AFTAC, Project  VT/007 
Final  Report,  Phase  1,  ARPA  Order  No.  104-50. 

D.  METHOD 


See  reference 


SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 

Title:  Theoretical  Isotropic  Processor 
COOP  Identification:  G640  ISOFIL 
Category:  Time  Series  Analysis 
Programmer:  D.  McCowan 

Date:  August,  1967 

B.  PURPOSE 


This  program  computes  and/or  applies  a  multichannel  iso¬ 
tropic  processor  to  seismic  array  data.  An  annular  ring 
noise  model  and  either  a  point  or  a  disc  signal  model  can  be 
specified.  The  program  then  solves  the  multichannel  Wiener- 
Hopf  equation  in  the  frequency  domain  to  get  the  optimum 
filter  which  rejects  the  noise  and  passes  the  signal.  The 
filter  is  written  on  a  save  tape  for  future  use.  An  option 
exists  to  filter  a  given  piece  of  data  and  plot  the  filter¬ 
ed  trace  and  the  direct  sum.  PERSON  plots  are  also  generated 
and  plotted  which,  in  conjunction  with  the  output  plots,  allow 
the  user  to  calculate  signal  to  noise  ratios. 

C.  USAGE 

1.  Operational  procedure:  This  is  a  FORTRAN- 6 3  Main 

program  with  the  following  subroutines:  DISMAT,  NMAT, 
SMAT1 ,  SMAT2 ,  RMAT1 ,  RMAT2 ,  APDEM,  INVERT,  MATMPY ,  GETEM, 
COOLCON,  PLOTEM,  BES ,  SPLOT,  DETRND,  and  PERSON66 .  In 
addition,  the  following  utility  subroutines  are  assum- 


ed 

to  be  on  the 

system  tape:  ERASE,  and  COOL. 

Parameters : 

a) 

NCH 

The  number  of  channels  in  the  filter, 
NCH<21. 

b) 

SR 

The  desired  sampling  rate,  usually  the 
same  as  the  data. 

c) 

SNR 

The  desired  signal  to  noise  ratio.  We 
have  found  SNR=4.0  to  work  well. 
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d )  VSG 


The  signal  velocity  for  a  disc  model 
if  one  is  used.  The  filters  will  pass 
"signals"  with  any  velocity  greater 
than  VSG. 

The  lower  velocity  limit  of  the  noise 
band . 


f )  VNH 


g)  IMODE 


h)  L 


i )  IPER 


j)  IDISCSIG 


k)  ISEIS 


The  upper  velocity  limit  of  the  noise 
band.  The  filters  will  reject  "noise" 
between  these  two  velocities. 

A  control  switch. 

./•  =  1  to  compute  the  filter  and  filter 
the  data 

=  0  to  read  the  filters  from  the  save 
tape  and  filter  the  data. 


The  number  of  points  in  the  frequency 
representation  of  the  filters.  This 
must  be  a  power  of  two  plus  one;  ie., 
17,  33,  129  etc.  L£l29.  The  number 
of  points  in  the  time  representation 
of  the  filter,  the  impulse  response, 
will  be  2*  (L-l)-.l, 

The  percentage  of  white  noise  to  add 
to  the  signal  model.  The  white  noise 
stabilizes  the  computation  and  assures 
good  convergence  in  the  filters.  We 
have  found  that  two  or  three  percent 
is  usually  adequate. 

A  control  switch  for  the  signal  model. 
=  0  for  an  infinite  velocity  signal 
model. 

=  1  for  a  disc  signal  model  V>VSG. 

The  seismogram  number  of  the  data  to 
be  filtered. 


1 )  NPTS 


The  number  of  points  to  be  filtered. 
NPTS<  3  500. 


m)  ISPT 


The  first  point  of  seismogram  ISEIS 
to  be  filtered. 


The  detrend  switch  for  the  data  to 
be  filtered. 

=  -  1  for  no  detrending . 

=  0  to  remove  the  mean. 

=  1  to  remove  the  mean  and  linear 

trend 
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o)  ICH(I) 

p)  X(I) 

q)  Y(I) 

r)  SF (I ) 


An  array  giving  the  positions  of  the 
NCH  data  channels  to  be  filtered  on 
the  input  tape.  ICH(I)  =  4  indicates 
the  Itfi  channel  is  the  fourth  on  the 
input  tape.  The  entries  in  ICH  must 
be  in  increasing  numeric  order. 

The  X  coordinate  of  the  I  data  chan¬ 
nel  . 

til 

The  Y  coordinate  of  the  I  data  chan¬ 
nel  . 

A  scale  or  demagnification  factor  for 
the  Ith  data  channel.  If  SF(I)  =  0.0 
or  blank  it  is  set  to  1.0. 


3.  Space  required: 


26K  locations. 


4.  Temporary  storage:  none 

5.  Alarms:  Various  alarms  for  illegal  data  are  printed  out 


6.  Error  returns:  none 

7.  Error  stops:  A  stop  occurs  after  each  alarm. 

ft  Tan*  mountings :  Tape  unit  1  is  the  save  tape  which  should 
*  be  a  scratclTMtape  when  the  filters  are  being  computed  and 
an  input  tape  when  they  are  being  read.  Tape  unit  i 
the  input  subset  tape  with  the  data  to  be  filtered.  Tape 
5  is  the  plot  tape  and  tapes  2  and  4  are  scratch  tapes. 

9.  input  and  output  formats:  Input  cards. 

Parameters  (a)  through  (j)  appear  on  the  first  card  as 
FORMAT  (15,  F10.2,  E15.3,  3F10.2,  «5)  Parameters  (k) 
through  (n)  appear  on  the  next  card  as  FORMAT  (110,  615). 
Parameters  (o)  through  (r)  appear  on  the  last  cards  as 
FORMAT  (3(13,  2F6.2,  F6.2)). 


Save  tape: 

The  save  tape  consists  of  one  file  of  binary  data  written 
in  two  logical  records.  The  first  record  is  a  127  word 
label  written  as  follows: 


LAB ( 1 ) 
LAB( 2 ) 
LAB ( 3 ) 
LAB ( 4 ) 
LAB( 5 ) 
LAB ( 6 ) 


NCH 

LAB( 7 )  =  L 

SR 

LAB (8'  =  IPER 

SNR 

LAB ( 9 )  =  IDISCSIG 

VSG 

LA3  (10— 10+NCH-l ) 

VNL 

LAB ( 31 — 31+NCH-l) 

VNH 
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The  second  record  is  a  5376  word  array  with  the  filters 
stored  close-packed  as  follows: 


10. 

11. 


Fi(t)  =  F(I ,t ) ,  ( (F(I , t ) ,  1=1,  NCH) ,  t=l , 2* (L-l)-l ) 
Selective  jumps:  none 


Timing:  A13  channel  129  point  filter  computed  and  applied 
to  3500  points  of  data  requires  about  24  minutes  of  CDC 
1604B  time. 


12. 

13. 


Accuracy:  Single  precision 


Cautions  to  user:  _  Do  not  exceed  ranges  of  NCH  and  L  as 


14. 

15. 


no  error  message  is  given  for  these  signs. 

Equipment  configuration:  Standard  COOP  with  RELOCOM. 


References :  Seismometer  Array  and  Data  Processing  System, 
Project  Veal-Uni form-AFTAC,  Project  i/T/007,  Final  Report,’ 
Phase  1,  ARPA  Order  No.  104-60. 


D.  METHOD 

See  reference. 
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SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 


Title:  Multichannel  Filter  Evaluation  Program 

COOP  Identification:  G639  LSTCHNCE 
Category:  Time  Series  Analysis 

Programmer:  D.  W.  McCowan 

Date:  November,  1967 

B.  PURPOSE 


This  program  computes  and  displays  the  f-k,  white  noise, 
and  amplitude  response  functions  for  theoretical  and  measured- 
noise  and  also  maximum  likelihood  filters.  The  programplots 
contours  of  the  f-k  transfer  function,  and  plots  the  white 
noise  and  amplitude  transfer  functions. 

C.  USAGE 


1.  Operational  Procedure:  This  is  a  Fortran-63  main  pro¬ 
gram  assuming  the  following  subroutines  on  the  system  tape: 
ERASE,  DISC63 ,  ABMAX ,  SKIPREC,  and  COOL. 


2.  Parameters : 


a)  JWORD 


b)  ISMPT 


c)  NF 


d)  VELR 


a  word  indicating  the  type  of  save 
tape  to  process. 

JWORD  =  6HMAXLIK  for  maximum- likeli¬ 
hood  filter 

JWORD  =  6HIS0FIL  for  theoretical  MCF 
JWORD  =  8HHEFALUMP  for  measured- 
noise  MCF 

the  number  of  points  in  frequency 
to  which  to  smooth  the  MCF  response 
functions.  This  must  be  a  power  of 
two  plus  one.  ISMPT<_6  5. 

the  number  of  frequencies  for  which 
f-k  contour  plots  of  the  f-k  trans¬ 
fer  function  are  desired  NF£l2. 

the  velocity  range  for  all  f-k  plots. 


e)  I LOG 


f)  NCODE 


-  a  switch  for-  linear/log  plots  of  the 
white  noise  and  amplitude  transfer 
functions . 

ILOG  =  0  generate  linear  plots 
ILOG  i  0  generate  db  plots 

-  the  spacing  between  plotted  points  in 
hundredths  of  an  inch 


g)  RANGE  -  the  range  in  inches  for  all  plots 

h)  X  CD  -  the  X  coordinate  of  the  Ith  seis¬ 

mometer  . 

i)  Y  -  the  Y  coordinate  of  the  Ith  seis¬ 

mometer. 

j)  FREQ  (I)  -  the  Ith  requested  frequency 

Space  required:  32  K 

Temporary  storage:  disc 

Alarms:  various  alarms  for  illegal  data  are  printed 
out . 


Error  returns :  none 


Error  stops :  a  stop  occurs  after  each  alarm 


Tajpe  mountings:  Tape  unit  1  is  a  scratch  tape.  Tape 
unit  2  is  the  input  save  tape  being  processed.  Tape 
unit  3  is  the  output  plot  tape. 


Input  and  output  formats:  Parameter  (a)  appears  on 
the  first  card  as  FORMAT  (A8).  Parameters  (b)  through 
(g)  appear  on  the  next  card  as  FORMAT  (2(215,  F10.2)). 
Parameters  (h)  and  (i)  appear  on  the  next  cards  as 
FORMAT  ( 3 ( 3X,  2F8.2,  6X)).  Parameter  (j)  appears  on 
the  last  card  as  FORMAT  (12F6.2). 


Selective  jumps:  none 

Timing:  same  as  FKSPTRUM  for  time  series  the  length 
of  the  filters. 


Accuracy:  single  precision 

Cautions  to  user:  none 

Equipment  configuration:  standard  COOP  with  RELOCOM, 
three  (3)  tape  units  and  818  disc  file. 

References:  D.  W.  McCowan,  Multichannel  Filter 
Report,  In  preparation.  - 
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SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 

Title:  Maximum  Likelihood  Filter  Program 
COOP  Identification:  UES  Z124  MAXLIK 


Category:  UES  General 
Programmer:  D.  W.  McCowan 

Date:  June,  1967 


B.  PURPOSE 

This  FORTRAN-63  program  computes  and/or  applies  a  21 
or  30  point  maximum-likelihood  filter  (realizable  or  sym¬ 
metric)  to  seismic  array  data.  Printed  output  is  produced 
as  well  as  a  plot  tape  displaying  the  input  and  output  data 
and  useful  computed  information.  In  addition,  a  save  tape 
is  generated  with  all  pertinent  data  about  the  operation 
written  in  an  easily  accessible  form.  A  provision  exists 
to  run  the  program  from  the  save  tape  processing  new  data 
with  both  kinds  of  maximum-likelihood  filters.  In  this  case, 
either  filter  is  efficiently  re-computed  from  information 
written  on  the  save  tape  and  a  new  save  tpae  is  generated 
with  details  of  the  new  filter  and  its  operation. 


The  printed  output  includes  the  actual  summation  constraint 
satisfied  by  the  filter  coefficients.  This  is,  for  realizable 
filters , 


L 


I 

i=l 


fi(t) 


=  6 


t 

1 


and  for  symmetric  filters. 


L 


I 

i=l 


f±(t> 


The  user  can  see  from  this  how  accurate  the  Levinson  recur¬ 
sion  was  in  calculating  the  filters.  Also  included  are  partial 
energy  sums  for  both  output  traces,  the  maximum-likelihood  and 
the  direct  sum. 
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They  are: 


t  0 

P(t)  =  I  x(t)2  , 
t=l 

and  are  referred  to  in  the  output  as  PERSON  values.  These 
are  useful  for  computing  RMS  noise  measurements  over  some 
arbitrary  interval.  The  two  maximum  absolute  values,  one 
for  the  input  traces  and  the  other  for  the  two  outputs ,  are 
also  printed. 

The  plot  tape  contains  plots  of  the  input  traces,  all  to 
the  same  scale  factor  and  the  two  output  traces  also  to  their 
own  common  scale  factor.  The  PERSON  values  are  plotted  to 
separate  scale  factors. 

The  user  is  referred  to  the  various  papers,  both  published 
and  unpublished,  on  the  theory  of  maximum-likelihood  filters 
for  a  thorough  technical  discussion.  These  references  also 
provide  definitions  for  the  quantities  described  in  this  write¬ 
up.  This  particular  write-up  will  not  concern  itself  with 
fhe  theory  of  maximum— likelihood  filters,  the  general  practices 
of  maximum- likelihood  filtering,  or  the  operation  of  the  sub¬ 
routines  in  this  program  set.  For  the  rest  of  these  suojects 
the  user  is  also  referred  to  the  papers  and  reports  listed 
as  well  as  the  write-ups  of  the  mathematical  subroutines  we 
have  written  for  this  program. 

C.  USAGE 


Operational  procedure:  This  is  a  FORTRAN- 6 3  main  pro¬ 
gram  using  the  overlay  option  provided  with  the  1604-B  COOP 
monitor.  It  is  usually  run  from  the  relocatable  binary  deck 
under  standard  monitor  control;  although  it  can  be  run  from 
the  overlay  tape  using  the  LOADMAIN  system. 

In  general  two  kinds  of  runs  can  be  made  with  this 
program.  The  first  is  a  '’Compute"  run  and  entails, the  com¬ 
putation  of  the  maximum— likelihood  filters  from  the  ground  up. 
The  other  kind  of  run  is  a  "Filter-Only"  run  in  which  the 
fflfer  is  read  from  the  save  tape  and  applied  to  new  data. 
Because  of  this,  the  Parameters  section  is  divided  into  two 
s ,  the  first  describing  the  parameters  necessary  for  a 
"Compute"  run  and  the  second  part  describing  those  required 
for  a  "Filter-Only"  run. 

2.1  Parameters  -  compute  run; 

a)  IWORD  -  filter  length  switch, 

IWORD  =  8HMAXLIK21,  for  a  21  point 
filter 
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IWORD  =  8HMAXLIK30,  for  a  30  point 
filter. 

b)  ISEIS  -  the  seismogram  number  of  the  data 

from  which  the  filter  is  to  be  com¬ 
puted  . 

c)  NCH  -  the  number  of  channels  to  be  used  in 

computing  the  filter,  (NCH<_21  for  a 
21  point  filter,  NCHO.O  for  a  30  point 
filter. 

d)  NPTS  -  the  number  of  points  in  the  fitting 

interval,  (NPTS<_6000) 

e)  ISPT  _  the  first  point  of  the  fitting  inter¬ 

val  on  seismogram  ISEIS. 

f)  IDT  -  the  trend  switch  for  the  fitting  inter¬ 

val  data, 

IDT  =  -1,  no  detrending 
IDT  =  0,mean  removed 

IDT  =  1,  mean  and  linear  trend  removed. 

g)  JUMP  _  the  number  of  digital  points  per  unit 

correlation.  This  determines  the  samp¬ 
ling  rate  of  the  correlation  functions, 
and  thereby  the  filters,  in  terms  of 
the  sampling  rate  of  the  input  data. 

JUMP  =  1,  the  sampling  rate  of  the 

filters  is  the  same  as  that  of 
the  data. 

JUMP  -  2,  the  filters  are  to  be  sampled 
every  other  data  points. 

JUMP  =  3,  the  filters  are  to  be  sampled 
every  third  data  point. 

etc . 

h)  ISOR  -  the  filter  designation  switch, 

I50R  =  1,  compute  symmetrical  filter 
ISOR  =  0,  compute  realizable  filter 

i)  IFOC  -  the  run  designation  switch. 

IFOC  =  1,  "Filter-Only"  run 
IFOC  =  0,  "Compute"  run. 


j  )  ICH 


Array  of  positions  of  tne  NCH  data 
channels  to  be  re-subset  from  the 
input  SUBSET  tape;  i.e.,  NCH  =  3, 
ICH  =  (2,4,5)  indicates  that  the 
second,  fourth,  and  fifth  channels 
on  the  input  tape  are  to  be  used. 
The  channels  must  be  listed  in  in¬ 
creasing  order. 


k) 

SF 

Array  of  demagnification  on  scale 
factors  for  the  NCH  selected  data 
channels . 

1) 

ISEIS1  - 

the  seismogram  number  for  the  data 
to  be  filtered. 

m) 

NCH1 

the  number  of  channels  to  be  filtered 
NCH1  must  equal  NCH. 

n) 

NPTS1  - 

the  number  of  points  of  seismogram 
ISEIS1  to  be  filtered  (NPTSK6000 ) . 

o) 

ISPT1  - 

the  first  point  of  seismogram  ISEIS1 
tc  be  filtered. 

P) 

IDT1 

the  detrend  switch  for  the  data  to  be 
filtered. 

q) 

ICH1 

array  of  channel  positions  of  the  in¬ 
put  tape  for  the  data  to  be  filtered. 

r) 

SF1 

array  of  demagnification  factors. 

Parameters  - 

Filter-only  run: 

a) 

IWORD  - 

the  filter  length  designation  switch 

b) 

ISOR 

the  filter  designation  switch 

c) 

IFOC 

the  run  designation  switch. 

d) 

ISEIS1  - 

as  above 

e) 

NCH1 

as  above 

f ) 

NPTS1  _ 

as  above 

g) 

ISPT1  _ 

as  above 

h) 

IDT1 

as  above 

i) 

ICH1 

as  above 

j) 

SF1 

as  above 

-32- 


32  K 


3 .  Space  required: 

4.  Temporary  space:  none 

5.  Alarms :  Various  alarms  caused  by  illegal  data  requests 
are  flagged  by  subroutine  MAKESUB.  The  user  is  referred  to 
the. write-up  of  subroutine  MAKESUB  for  their  meanings.  In  ad¬ 
dition,  as  subroutine  MAKESUB  does  not  stop,  the  main  program 
also  indicates  an  error  has  taken  place.  It  is  flagged,  though 
not  designated,  and  the  job  is  terminated. 

6.  Error  returns :  none 

7 •  Stops :  The  program  halts  if  an  error  condition  in 
subroutine  MAKESUB  is  detected. 

8.  Tape  mountings :  The  input  subset  tape  for  either  kind 
of  run  must  be  mounted  on  logical  tape  1.  The  plot  tape  is 
logical  tape  4  and  the  save  tape  is  logical  tape  3.  In  a  "filter 
only"  run,  another  save  tape  must  be  input  on  logical  tape  5. 
Logical  tape  5  is  used  only  for  scratch  in  a  "compute"  run. 
Logical  tape  6  is  the  overlay  tape  and  must  be  low  density  (the 
cover  card  must  contain  this  instruction).  Logical  tape  2  is 
a  scratch  and  can  be  equivalenced  to  56,  as  it  is  initially 
rewound. 


9*1  Card  formats  -  compute  run:  Parameter  (a)  appears 
on  the  first  card  as  FORMAT  (A8).  Parameters  (b)  through  (i) 
appear  on  the  second  card  as  FORMAT  (8110).  Parameter  ( j ) 
appears  on  the  next  card  as  FORMAT  (2113).  Parameter  (k)  ap¬ 
pears  on  the  next  card  as  FORMAT  (4E20.3).  Parameters  (1) 
through  (p)  appear  on  the  next  card  as  FORMAT  (5110).  If 
parameter  (1)  is  zero  or  blank,  no  further  cards  are  read  and 
the  data  to  be  filtered  will  be  the  same  as  the  data  used  to 
compute  the  filter.  Otherwise,  parameter  (q)  appears  on  the 
next  card  as  FORMAT  (2513)  and  parameter  (m)  appears  on  the 
last  cards  as  FORMAT  (4E20.3). 

9.2  Card  formats  -  filter-only  run:  Parameter  (a)  appears 
on  the  first  card  as  FORMAT  (A8).  Parameters  (b)  and  (c)  ap¬ 
pear  on  the  second  card  as  FORMAT  (60X,  2110).  Parameters  (d) 
through  (h)  appears  on  the  next  card  as  FORMAT  (2513)  and  para¬ 
meter  (i)  appears  on  the  next  card  as  FORMAT  (5110).  Parameter 
( j )  appears  on  the  last  cards  as  FORMAT  (4E20.3). 

9 • 3  Save  tape  format:  The  save  tape  consists  of  8  files 
of  binary  data  as  described  below: 

File  Record  Length  Quantity  Description 

1  1  225  label 
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LAB  (1)  =  ISEIS 
LAB  (2)  =  NCH 


File  Record  Length 


Quantity 


EOF 


Description 


LAB 

(3) 

a 

NPTS 

LAB 

(4) 

a 

ISPT 

LAB 

(5) 

a 

IDT 

LAB 

(6) 

= 

JUMP 

LAB 

(7) 

a 

ISOR 

LAB 

(8) 

a 

IFOC 

(LAB 

(39+1) ,1=1 

NCH)=ICH 
(LAB  (39+1), 1=1, 
NCH)=SF 
LAB  (100)  =  L 


2 

EOF 

1 

9261 

3 

1 

255 

EOF 

4 

EOF 

1 

9261 

5 

EOF 

1 

441 

R 

label 

E=R'1 IIG(GTR-1G)_1 

F=ED 


See  Reference 

Same  as  above  but 
for  data  to  be 
filtered . 

See  reference  3. 

See  reference  3. 


EOF 

7 


I  127  or  256  label 

2-NCH+l  LAB  (3)  data 

II  LX 


EOF 

8 

EOF 


1 

2 


LX 

1 

LX 


Phased  sum 

LX 

max-like 


Input  tape  label 
as  set  for  data  to 
be  filtered,  data 
channels 

Length  of  phased 
sum 


Length  of  max- like 
trace 


10.  Selective  jumps:  none 

11.  Timing:  The  times  required  to  calculate  the  correlation 
matrix  and  the  filter  are  given  by: 

T  =  1.83  x  10"6.N2.L.LX 
c 

Tf  =  1.65  x  10"5.N3.L2 
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respectively.  T  is  in  minutes  of  1604-B  time,  N  is  the  number 
of  channels,  L  is  the  filter  length,  and  LX  is  the  number  of 
points  in  the  fitting  interval.  The  time  +o  do  the  convolut¬ 
ions  and  generate  the  plot  varies;  however,  19  channels  of 
4800  points  each  requires  about  12  minutes  of  1604-B  time. 

12.  Accuracy :  Single  precision. 

13.  Caution  to  users:  Text  above 

14.  Equipment  configuration:  Standard  COOP.  Typical 
COOP  cards  for  a  "Compute1'  run  and  a  "Finite-only"  run  are 
given  below: 

a  "Compute"  run  and  a  "Filter-Only"  run  are  given  below: 

7  COOP,  23499,  BULOVA,  1/ 1/S/ 3/4/ 5/ 6/ 56/57/E/ 2  =  56,  140, 

9  500. 


7  COOP,  23499,  BULOVA,  1/1/ 5/S/ 3/4/ 5/ 6/56/57/E/ 2  =  56, 

9  140,  500. 

15 .  References : 

(1)  Bruffey,  V.,  Write-up  of  Z122  MAKESUB,  Seismic 
Data  Laboratory,  March,  1967. 

(2)  Fletcher,  N.H.  Write-Up  of  Z24  SUBSET,  Seismic 
Data  Laboratory,  September,  1964. 

(3)  Flinn,  E.  A.,  and  Claerbout,  J.F.,  Some  Topics 
in  Digital  Filtering,  Theory  and  Application, 
Seismic  Data  Laboratory,  Internal  Memorandum, 
September  16,  1965. 

(4)  Flinn,  E.  A.,  et.al.,  Two  Examples  of  Maximum 
Likelihood  Filtering  of  LASA  Seismograms, 
Seismic  Data  Laboratory,  Report  No.  148,  June, 
1966  . 

(5)  Flinn,  E.A.,  et.al.,  Maximum  Likelihood  Filter¬ 
ing  of  LASA  Noise  Seismograms,  Seismic  Data 
Laboratory,  Report  No.  149,  June  8,  1966. 

(6)  McCowan,  D.  W. ,  Write-Up  of  E509  PHILTRE, 
Seismic  Data  Laboratory,  November,  1965. 

(7)  McCowan,  D.  W.,  Write-Up  of  Z78  MCURT,  Seismic 
Data  Laboratory,  September,  1965. 

(8)  McCowan,  D.W.,  Write-Up  of  E509  MAX-LIKE, 
Seismic  Data  Laboratory , March,  1965. 

D.  METHOD 


See  references. 
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SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 

Title:  Multi-channel  Filter  to  Paper-Tape  Format 

COOP  Identification:  M220  MCFONPT 

Category:  Conversion 

Programmer:  C.  W.  Clowes 

Date:  2  October  1967 

B.  PURPOSE 

Reads  either  on  ISOFIL,  MAXLIK  or  HEFALUMP  save  tape. 
Scales  and  converts  selected  channels  to  punch  paper  tape 
format  and  places  on  a  MCFONPT  save  tape.  The  MCFONPT  save 
tape  is  punched  by  16 OA  program  PUNCH  MCFONPT. 

C.  USAGE 


1*  Operational  Procedures :  Operator  must  be  supplied 
with  3  data  cards  describing  tape  to  read,  channels  on  tape 
and  number  of  filters  to  punch. 

2 .  Parameters : 

a)  Card  Format  (13,  6X,  A8,  2X,  2I3/(BI3)) 

b)  Card  1 

cc  1-3  NF  -  Number  of  filter  channels  to  be  punched. 
10-17  ISW  -  Name  of  input  tape,  (left  justified) 
MAXLIK,  ISOFIL  or  HEFALUMP 
20-22  SC  -  Scaling  constant  (usually  63) 

23-25  IOIND  -  Output  channel  index. 

c)  cc  1-39  Array  ICH(I),  1=1,13 

Array  for  channel  identifiers  on  input  tape. 
Range  0-12.  Puncn  starting  column  1,  3 
columns  per  channel  number,  is  ascending 
order.  Example:  ICH(?)=4,  Indicates  2nd 
filter  on  magnetic  tape  corresponds  to  the 
4th  array  channel,  i.e.  Z4. 

d)  cc  1-39  Array  JCH(I),  1=1, B. 

Array  identifying  channels  that  are  to  be  punched. 
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Range  0-12.  Punch  on  card  starting  column  1,  3 
columns  per  channel  number,  in  ascending  order. 
Caution:  Channels  listed  in  this  array  must  be 

listed  in  ICH  array. 

3.  Space  Required:  Program  requires  14600  locations 

4 .  Temporary  Storage: 

5 .  Print-Outs : 

a)  Channel  No.  selected  is  not  on  input  tape  or 
data  card  is  wrong . 

b)  Printout  also  lists: 


i . 

input  tape  name 

ii . 

number  of  filters 

on  tape 

iii . 

number  of  filters 

to  be  processed 

iv . 

channel  numbers  on 

tape  (ICH) 

v . 

channel  numbers  to 

process  (JCH) 

c)  A  compiled  list  of  all  of  the  filters  processed 
is  printed. 

S.  Error  Return:  If  an  error  in  the  number  of  filter 
channels  selected  occurs,  control  is  returned  to  the 
monitor. 

7.  Error  Stops :  None 

8.  Input  and  Output  Tape  Mountings:  Input:  Lo  Unit  1 

. .  Output:  Lo  Unit  2 

9 .  Input  and  Output  Formats: 

a)  Card  input  (see  parameters) 

b)  Tape  input  (1)  is  output  from  MAXLIK  (Unit  3), 

ISOFIL  (Unit  1)  HEFALUMP  (Unit  3). 

c)  Tape  output  (2).  13  records  of  13  channels.  Each 

record  consists  of  one  channel  x  number  of  filter 
coefficients . 

10.  Selective  jump  and  stop  settings:  Switch  1  on  to 
bypass  printout  of  filter  coefficients . 

11.  Timing:  Processing  maximum  of  13  channels  of  255 
filters  coefficients  with  printout  takes  6  minutes. 

SW  1  on  takes  considerably  less. 

12.  Accuracy:  Filters  are  scaled  and  truncated  to  a  four 
digit  octal  integer,  range  -3777  to  +3777. 

13.  Caution  to  User:  Channel  filters  range  from  0  through 
12' on  data  cards  and  printouts.  Filters  listed  on  card 
3  must  be  present  on  tape  and  card  2. 
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14 •  Equipment  Configuration:  Standard  COOP.  1604  pro¬ 
cessor,  2  tape  drivers,  1607  printer,  card  reader. 

15 •  References :  Write-up  of  HEFALUMP,  ISOFIL,  MAXLIK 

by  D.  W.  McCowan.  Texas  Instruments,  1966,  Operation 
and  Maintanence  of  The  CPO  Multichannel  Filter  System. 
Instruction  Manual. 

D.  METHOD 

The  data  is  read  from  either  tapes  ISOFIL,  MAXLIK  or 
HEFALUMP  into  array  FILCOEF.  Channels  to  be  processed  (JCH) 
are  compared  to  channels  on  tape  (ICH)  and  the  channels  listed 
in  JCH  are  scaled.  All  scaled  filters  are  then  compared  using 
the  following  scheme: 

NF  •  LF  •  SC  •  FMAX  =  A 


where 


NF  is  number  of  filter  channels  to  punch. 
LF  is  number  of  filter  points/channel . 

SC  is  scale  constant  (from  data  card) 
FMAX  is  max.  absolute  filter  point. 


if  A  >  16.000,0001Q,  all  filter  points  are  scaled  again  by 

16  X  106/A. 


All  channels  (0-12)  are  placed  in  160A  punch  tape  format.  For 
those _ channels  not  listed  in  JCH  (punch  card),  the  filter  co¬ 


efficients  are  set  at  zero, 
format  of: 


Each  filter  coefficient  takes  the 


Word  1 

Word  2 

Word  3 

12  !  IT” 

i  bits  :  bits 

mi 

12  |  TT 

l 

!  12  m  i2~y 
!  !  !  1 

« _ _ 

V _  _ t 

Channel 
index  (BCD) 

Filter  coefficient 
index  (BCD) 

mem. 

sub 

output 

channel 

Filter 

Coefficient 

space 

(always 

1) 


index 
(from 
card ) 


(octal) 


AH  of  the  filter  coefficients  for  eacn  channel  are  combined 
into  one  record,  and  placed  on  magnetic  tape.  There  are  13 
records  on  the  output  tape  ( 2 ) . 


SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 

Title:  Multi-channel  Filters  to  Paper  Tape 

SWAP  Identification:  M221  Punch  MCFONPT 
Category:  Conversion 

Programmer:  C.  W.  Clowes 

Date:  2  October  1967 

B.  PURPOSE 

This  is  a  160-A  program  using  1  tape  unit  and  tape  punch. 
Reads  save  tape  from  program  MCFONPT  and  punches  same  on  8 
channel  paper  tape.  Must  use  1  inch,  eight  channel  tape,  with 
7  punches. 

C.  USAGE 

1.  Operational  Procedure:  Load  Bi-octal  tape  into  blank 
0.  Turn  on  punch,  (check  tape  supply),  master  clear  and  run. 
Program  will  read  input  tape  (1)  and  punch  onto  paper  tape. 
When  finished,  tape  will  re-wind  and  stop  at  P=0214,  A=0066. 
Run  is  completed.  To  re-run,  master  clear,  run. 

3.  Space  required:  Program  uses  113  160A  words. 

7.  Error  stops:  P  =  0140  A  =  0077.  Unable  to  read  chan¬ 

nel  record,  set  jumpt  key  1  to  try  reading  again,  off, 
process  with  error. 

8.  Input  and  Output  Tape  Mountings:  Unit  1.  Input  tape. 

9.  Input  and  Output  Format:  Unit  1  input,  13  records  of 

filter  coefficients.  1  record  per  channel. 

10.  Selective  jump  switches:  (see  7,  error  stops) 

11.  Timing:  About  3  minutes 

13.  Caution  to  User:  Make  sure  1"  tape  is  in  punch. 

14.  Equipment  Configuration:  160  A  computer  with- paper 
tape  punch  and  CDC  162  magnetic  tape  controller. 
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METHOD 


3  feet  of  leader  (177g)  is  placed  on  paper  tape  followed 

by  13  records,  read  from  magnetic  tape  (1)  and  punched  one  at 
a  time.  Each  record  is  separated  on  paper  tape  by  1  tape 
feed  character.  3  feet  of  trailer  is  punched  on  the  end  of 
paper  tape . 


SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 
Title:  MULTICOH 

COOP  Identification:  G630  MULTICOH 
Category:  Time  Series  Analysis 

Programmer:  J.  Jih 

Date:  August,  1967 

B.  PURPOSE 

This  program  computes  multiple  coherence  functions  for 
seismic  array  data  rapidly  and  efficiently.  Given. an  ori¬ 
ginal  set  of  N  subset  data  channels,  the  program  will  com¬ 
pute  the  N-l  multiple  coherence  functions: 

Xi(N-i/N, - ,N-i+l)  i  =  1...N-1 

The  program  will  then  reorder  the  N  data  channels  any  number 
of  times,  each  time  computing  another  N-l  multiple  coherence, 
function.  The  print-out  includes  a  description  of  the  notation 
used.  Optional  print-out  includes  all  the  auto. and  cross  spectra. 
In  addition  a  provision  exists  to  plot  the  multiple  coherence 
functions.  They  Cooley-Tukey  method  of  spectral  estimation 
is  used  to  botain  high  speed, 

C.  USAGE 

1.  Operational  procedure :  This  is  a  Fortran-63  main  pro¬ 
gram  with  the  following  subroutines:  SUBTAPE,  DETRND,  XMNG1, 
COOLER,  SMOOTH,  PRINTSP ,  NOTATION,  RESCALAT,  CP LOT,  and  SPLOT . 

In  addition,  the  following  utility  subroutines  are  assumed  to 
be  on  the  system  tape:  SKIPREC,  ERASE,  DISC63,  and  COOL. 

2.  Parameters:  This  program  is  set  up  to  do  multiple 
runs.  It  reads  cards  until  an  EOF  is  encountered  and  then 
terminates . 

a)  ISM 

b)  N 


The  seismogram  number  of  the  desired  data. 

The  number  of  data  channels  to  be  used  from 
seismogram  ISM.  (2  <_  N  £  20) 
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c)  LX 


d)  L 


e)  ISPT 

f)  IDT 


g)  INDX 


h)  MP 


i)  IPLOT 


The  number  of  digital  points  to  be  used 
from  each  data  channel.  LX  must  be  a 
power  of  two  (LX  <  4096).  If  LX  >  4096 
it  will  be  truncated  to  4096. 

The  number  of  points  desired  in  the  mul¬ 
tiple  coherence  functions  between  DC  and 
the  folding  frequency.  L  must  be  a  power 
of  two  plus  one,  i.e.,  17,  33,  65,  129 
etc.  L  must  also  satisfy  the  following 
bounds. 

3  £  L  £  33  for  15  <  N  <_  20 

3  £  L  £  65  for  11  <  N  <  15 

3  £  L  £  129  for  7  <  N  _£  11 

3  £  L  £  257  for  5  <  N  <  7 

3  £  L  £  513  for  2  £_  N  £  5 . 

if  these  conditions  are  not  met,  L  will 
be  truncated. 

The  first  digital  point  of  the  desired 
data.  (ISPT  £  1) 

The  detrend  switch: 

0  or  blank  to  remove  the  mean 

1  to  remove  the  mean  and  linear 

trend 


for  no  detrending 

The  number  of  times  to  change  the  order 
of  the  original  data  channels  and  re¬ 
compute  the  coherence  functions.  In 
all,  INDXMN-l)  multiple  coherence  func¬ 
tions  will  be  computed  because  the  ori¬ 
ginal  order  will  not  be  changed. 

Switch  for  printing  the  power  spectra 

=  0  or  blank  for  no  printout 


*  "to  print  power  spectra 

Switch  for  plotting  the  multiple  co¬ 
herence  functions 
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0  or  blank  for  no  plot 


j) 


k) 

l) 


m) 


3. 


4. 


5. 


6. 


7. 


8. 


9. 


10. 


11. 


t  0  to  plot  coherences. 

The  spacing  between  plotted  points  in 
hundredths  of  an  inch.  All  plots  are 
to  a  five  inch  range. 

The  position  of  the  Ith  data  channel  on 
the  input  tape. 

The  demagnification  factor  of  the  Ith 
data  channel.  SF(I)  is  divided  into 
every  point  of  the  Ith  channel.  If 
SF(I)  is  0.0  or  blank  it  is  set  to  1.0. 

An  array  of  channel  positions  giving  a 
new  order  for  the  data  channels  and  hence 
new  combinations  of  multiple  coherencies. 
The  N  numbers  in  this  array  must  be  the 
same  as  those  in  ICH  but  can  be  in  any 
order. 

Space  requ'red:  20190  locations. 

Temporary  storage:  DISC 

Alarms :  Various  error  messages  are  printed  out  for 

errors  in  the  input  data  and  parameters. 

Error  returns:  none 

Error  stops:  There  are  no  stops.  The  program  always 
proceeds  with  the  next  job  request. 

Tape  mountings ;  An  input  tape  of  SUBSET  seismograms 
must  be  on  logical  unit  1.  A  scratch  tape  must  be 
mounted  on  lcgical  unit  2.  If  plots  are  requested, 
an  output  tape  must  be  mounted  on  logical  unit  3. 

Formats :  For  each  case  to  be  run  the  card  deck  is 

as  follows:  The  first  card  lists  parameters  (a) 
through  (j)  as  FORMAT  (1018).  The  second  card  lists 
parameters  (k)  and  (1)  as  FORMAT  (5(12, F14. 2) ) .  If 
INDX=0  no  further  cards  are  read  for  this  case.  If 
INDX?0  the  next  INDX  cards  list  parameter  (m)  as 
FORMAT  (2013).  The  deck  for  the  next  request  immed¬ 
iately  follows. 

Selective  jumps:  none 

Timing:  A  N=5,  LX=i024,  L=17,  INDX=0,  MP=1,IPL0T=1, 
case  takes  4  minutes  of  CDC  1604B  time. 


NCODE 

ICH(I) 
SF  ( I ) 

IC(I) 
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12.  Accuracy :  Single  precision 

13.  Cautions  to  user:  none 

l4-  Equipment  configuration:  Standard  COOP 
15.  References: 

a)  Bendat ,  and  Piersol,  Measurement  and  Analysis 
of  Random  Data,  John  Wiley  S  Sons,  New  York,  1966. 

b )  Frazer,  Duncan,  and  Collar,  Elementary  Matrices, 
Cambridge  University  Press,  Cambridge,  1963,  pp. 11 2-11 8. 

c)  Jih,  write-up  of  subroutine  XMNG ,  SDL,  1966. 

D.  METHOD 


The  multiple  coherences  of  channel  1  with  channels  2  and 
3  removed  is  computed  from  the  augmented  spectral  matrix  S': 


i 


where  S  is  the  spectral  matrix  of  the  channels  to  be  removed 
and  G  is  the  vector  of  cross-spectra  between  them  and  the  de¬ 
sired  channel.  The  formula  is: 

A2  ( 1/  2  .  3  )  =  1  -  '1 - 

Sii  det  S 

The  ratio  of  the  determinants  is  calculated  directly  from  the 
escalator  method  of  matrix  inversion.  It  is  easily  seen  that 
all  of  the  multiple  coherences : 

Ai(N-i/N, . . . ,N-i  +1)  i  =  1  ...  N-l 

can  be  calculated  by  the  inversion  of  S'  by  saving  the  sub¬ 
ratios  of  determinants  in  the  escalator  process. 
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SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 
Title:  PARTLCOH 

COOP  Identification:  G631  PARTLCOH 
Category:  Time  Series  Analysis 

Programmer:  J.  Jih 

Date:  August,  1967 

B.  PURPOSE 

This  program  computs  partial  coherence  functions  for  taped 
data.  The  program  also  computes  the  amplitude  and  phase  of 
the  associated  transfer  function.  The  output  includes  print¬ 
outs  of  these  functions  as  well  as  plots  of  the  coherence 
functions.  The  Cooley-Tukey  method  of  spectral  estimation 
is  used  to  obtain  high  speed. 

C.  USAGE 

1.  Operational  procedure:  This  is  a  Fortran-63  main  pro¬ 
gram  with  the  following  subroutines:  SUBTAPE,  DETRND,  XMNG2 , 
COOLER,  SMOOTH,  BP LOT,  PSP,  RESCALAT ,  and  SPLOT.  In  addition 
the  following  utility  subroutines  are  assumed  to  be  on  the 
system  tape:  SKIPREC,  ERASE,  DISC63,  and  COOL. 

2.  Parameters:  This  program  is  set  up  to  do  multiple 
runs.  It  reads  cards  until  an  EOF  is  encountered  and  then 
terminates . 

a)  ISM  -  The  seismogram  number  of  the  desired 

data 

b)  n  -  The  number  of  data  channels  to  be  used 

from  seismogram,  ISM.  (2  <  N  £  20). 

Any  channel  used  to  specify  a  partial 
coherence  function  must  be  a  member  of 
this  original  set. 

c)  LX  -  The  number  of  digital  points  to  be  used 

from  each  data  channel .  LX  must  be  a 
power  of  two  (LX  <_  4096).  If  LX  >  4096, 
it  will  be  truncated  to  4096. 
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d)  L 


The  number  of  points  desired  in  the 
partial  coherence  functions  between  DC 
and  the  folding  frequency.  L  must  be 
a  power  of  two  plus  one,  i.e.,  17,33, 

65,  129,  etc.  L  must  also  satisfy  the 
following  bounds. 

3  <_  L  £  33  for  15  <  N  <  20 

3  <  L  <  65  for  11  <  N  <  15 

3  <_  L  <_  129  for  7  <  N  <  11 

3  <  L  <_  257  for  5  <  N  <  7 

3  <_  L  £  513  for  2  <  N  <_  5 . 

If  these  conditions  are  not  met,  L  will 
be  truncated. 

e)  ISPT  -  The  first  digital  point  of  the  desired 

data.  (ISPT  >  1) . 

f)  IDT  -  The  detrend  switch. 

s  0  or  blank  to  remove  the  mean 
=1  to  remove  the  mean  and  linear 

trend. 

=  -1  for  no  detrending 

g)  INDX  -  The  number  of  partial  coherence  functions 

to  compute.  (INDX  >_  1). 

h)  MP  -  Switch  for  printing  the  power  spectra 

=  0  or  blank  for  no  print-out 
^  0  to  print  power  spectra. 

i)  IPLOT  -  Switch  for  plotting  the  partial  co¬ 

herence  functions. 

=  0  or  blank  for  no  plot 
i  to  plot  coherences 

j )  NCODE  -  The  spacing  between  plotted  points  in 

hundredths  of  an  inch.  All  plots  are 
to  a  five  inch  range. 

The  position  of  the  Ith  data  channel  on 
the  input  tapes.  (I  -  1,2,...,N).  This 
specifies  the  original  set  of  channels. 


k)  ICH(I) 


1) 


I 


l)  SF(I)  -  The  demagnif  icati<  factor  for  the  Ith 

data  channel.  SFl  )  is  divided  into 
every  point  of  the  Ith  channel.  If 
SF(I)  is  0.0  or  blank  it  is  set  to  1.0. 

m)  NCH  -  The  number  of  data  channels  involved  in 

this  particular  partial  coherence  func¬ 
tion.  NCH  <  N. 

n)  IC(I)  -  An  "array  giving  the  channel  positions 

of  the  data  channels  that  specify  this 
particular  partial  coherence  function. 

The  first  two  numbers  specify  the  two 
channels  between  which  the  coherence 
is  to  be  calculated.  The  next  NCH-2 
numbers  specify  the  channels  to  be  pre¬ 
diction  filtered  out.  The  numbers  in 
this  array  must  be  the  same  as  those  in 
the  ICH  array  but  can  occur  in  any  order, 
numbt  or  combination.  (1  =  1,2, ...  ,NCH) . 

3.  Space  required:  20190  locations. 

4.  Temporary  storage:  DISC 

5.  Alarms :  Various  error  messages  are  printed  out  for 
jrrors  in  the  input  data  and  parameters. 

6.  Error  returns;  none 

7.  Error  stops;  There  are  no  stops.  The  programs  always 
proceeds  with  the  next  job  request. 

8.  Tape  mountings ;  An  input  tape  of  SUBSET  seismograms 
must  be  on  logical  unit  1.  A  scratch  tape  must  be 
mounted  on  logical  unit  2.  If  plots  are  requested, 
an  output  tape  must  be  mounted  on  logical  unit  3. 

9.  Formats :  For  each  case  to  be  run  the  card  deck  is 
as  follows: 


The  first  card  list  a  parameters  (a)  through  (j)  as 
FORMAT  (1018).  The  second  card  lists  parameters  (k) 
and  (1)  as  FORMAT  (5(12,  F14 . 3 ) ) .  The  next  INDX  cards 
list  parameters,  (m)  and  (n)  as  FORMAT  (2113).  The 
deck  for  the  next  request  immediately  follows. 

10.  Selective  jumps:  none 

11.  Timing:  A  N=5,  LX=1024,  L=33,  INDX=3,  MP=0,  IPL0T=1, 
NCH=5,4, 3,  case  takes  4  minutes  of  CDC  1604B  time. 
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12.  Accuracy :  Single  precision 

13.  Cautions  to  user:  none 

14 *  Equipment  configuration:  Standard  COOP 
1 5 .  References : 

a)  Bendat,  and  Piersol,  Measurement  end  Analys 
of  Random  Data,  John  Wiley  £  Sons,  New  York 

b)  Jih,  write-up  of  subroutine  XMNG,  SDL,  1966 
D  METHOD 

See  reference  (a),  Page  315,316. 


i_s 

1966  . 
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SEISMIC  DATA  LABORATORY 
Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 

Title:  F-K  Plotter 

COOP  Identification:  G636  PLOTFK 
Ca.egory:  Time  Series  Analysis 
Programmer:  D.  W.  McCowan 

Date:  October  1967 

B.  PURPOSE 

This  program  takes  a  spectral  matrix  read  from  a  binary 
program  HEFALUMP  save  tape  and  plots  requested  f-k  spectra. 
It  was  primarily  designed  to  plot  f-k  spectra  for  averaged 
spectral  matrices  computed  by  program  SPECAVE  but  will  work 
with  any  HEFALUMP  save  tape. 


C.  USAGE 
1.  0 


Operational  Procedure:  This  is  a  Fortran- 6 3  main  pro¬ 
gram  with  the  following  subroutines:  FKMAT,  CNTRU5, 
LOCATION,  READITIN,  and  SETCONST.  In  addition,  the 
following  utility  subroutines  are  assumed  to  be  on 
the  system  tape:  DISC63. 

Parameters : 


a) 

NF 

The 

number  of  frequencies  at 

which  f-k 

plots  are  desired.  NF  £  12. 

b) 

VELR 

The 

velocity  range  for  all  f- 

-k  plots. 

c) 

FREQ (I)- 

The 

Ith  requested  frequency. 

Space  Required:  32K  locations 
Temporary  Storage:  DISC 
Alarms:  Alarms  for  reque 


Alarms  for  requested  frequencies.  LE, 
the  DC  frequency  and  .GE.  the  folding 
frequency  are  printed. 


Error  Returns:  none 
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* 


* 


7 .  Error  Stops :  none 

8*  Tape  Mountings :  An  input  save  tape  generated  by 
either  program  SPECAVE  or  program  HEFALUMP  must 
be  on  unit  1. 

9.  Input  Format;  Parameters  (a)  and  (b)  appear  on  the 
first  card  as  FORMAT  (15,  F10.2).  Parameter  (c) 
appears  on  the  next  card  as  FORMAT  (12F6.2). 

10.  Selective  Jumps:  none 

11.  Timing :  A  NCH=13,  L=65,  NF=6  case  takes  16  minutes 
of  CDC  1604B  time. 

12.  Accuracy :  Single  precision 

13.  Cautions  to  User:  none 

l1*.  Equipment  Configuration:  Standard  COOP  with  RELOCOM 
15.  References ; 

Write-ups  for  HEFALUMP  and  SPECAVE  by  D.  W.  McCowan 
and  FKSPTRUM  by  J.  Jih. 

D.  METHOD 

See  References 


f 
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SEISMIC  DATA  LABORATORY 
Alexandria, Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 
Title:  PREDICT 

COOP  Identification:  G628  PREDICT 

Category:  Time  Series  Analysis 
Programmer:  V.  R.  Bruffey 

Date:  April  28,  1967 

B.  PURPOSE 

To  compute  a  prediction  of  the  straight  sum  of  multi¬ 
channel  arrays,  ispan  units  ahead  in  time  where  ispan  is  the 
prediction  span.  The  predicted  sum  and  the  actual  sum  are 
subtracted  to  yield  a  prediction  error,  and  all  are  plotted 
to  the  same  scale  factor.  Also,  to  compute  and  plot  power- 
spectra  of  these  traces. 

C.  USAGE 

1.  Operational  procedure:  The  program  is  written  in 

Fortran- 6 3  and  uses  the  following  subroutines: 

(a.)  MAKESUB  -this  subroutine  reads  an  input  seis¬ 
mogram  and  produces  another  seismo¬ 
gram  with  the  desired  data.  It  uses 
the  subroutine  DETRND  which  removes 
the  mean  or  mean  and  linear  trend. 

(b.)  MCORT  -computes  a  correlation  matrix.  It 

uses  subroutine  DOT  which  performs  a 
vector  product. 

(c.)  BLMKGC  -computes  a  prediction  filter  by  the 
Levinson- Robinson  recursion  method. 

It  uses  subroutines  MAINE  and  MATMPY 
which  perform  matrix  inversion  and 
matrix  multiplication,  respectively. 


(d.)  PHILTRE  -convolves  a  filter  with  data  to  ob¬ 
tain  filtered  data. 


(e.)  SPECKS  -  computes  the  power-spectra  of  an  array. 

It  uses  the  subroutines  SMOOTH  and  COOLER 
for  smoothing  data  and  computing  the 
Fourier  transform. 

(f.)  XPLOT  -  writes  a  plot  tape  with  a  given  scale 

factor . 


2.  Options :  The  program  can  be  run  three  ways,  a  Stand¬ 
ard  run,  an  Option  1,  and  an  Option  2.  The  list  below 
shows  the  calculations  performed  by  each  method. 


Standard  Method 
Computes  correlation  matrix 

Computes  prediction  filter 
Filters  data 
Writes  save  tape 
Plots  results 
Computes  power-spectra 
Plots  power-spectra 


Option  1 

Reads  correlation  matrix 
from  a  save  tape 

Computes  prediction  filter 

Filters  data 

Writes  save  tape 

Plots  results 

Computes  power-spectra 

Plots  power-spectra 


Option  2 

Reads  filter  from  save  tape 
Filters  data 
Plots  results 


Computes  power-spectra 
Plots  power-spectra 

Note  that  Options  1  and  2  assume  an  input  save  tape  generated 
by  an  earlier  run  of  the  program,  and  that  Option  2  does  not 
generate  a  save  tape.  A  description  of  the  save  tape  is  given 
under  the  section  Input  and  output  tape  mountings. 


3 .  Parameters : 

a)  ISEIS  - 

b)  NCH 

c)  NPTS 


the  seismogram  number  of  the  data 
the  number  of  channels  to  be  processed 
the  number  of  data  points  in  the  sample 
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d)  ISPT 


the  first  point  -'f  the  requested  data 
sample  for  each  channel  on  the  input 
tape . 


e) 


f ) 


g) 

h) 


i) 


j) 


k) 


1) 


m) 


n) 


4. 


IDETREND 

JUMP 


switch  for  detrending  the  data.  IDT=  -1 
for  no  detrending,  IDT  =  0  to  remove  mean 
and  IDT  ~  1  to  remove  and  linear  trend. 

the  number  of  elements  skipped  for  each 
lag  in  correlation  used  in  computing  the 
filter. 


I  SPAN 
LENGTH 


I0PTI0N 


I  CPI 


SF 


I START 


N 


NSM 


prediction  span  in  points 

the  number  of  lags  of  correlation  used  in 
computing  the  filter.  Length  is  the  num¬ 
ber  of  elements,  for  each  channel,  in 
the  filter. 

switch  for  selecting  the  option.  I0PTT0N 
0  or  blank  means  standard  run,  I0PTI0N  = 

1  or  2  means  Option  1  or  Option  2,  res¬ 
pectively. 

an  array  of  length  NCH  which  specifies 
the  data  channels  indices.  For  example, 
if  the  input  tape  contains  channels  Zl, 

Z2,  Z6,  Z7,  and  Z9,  respectively,  and 
one  wants  to  process  Zl,  Z6  and  Z9,  then 
ICH  (1)  =  1,  ICH  (2)  =  3,  and  ICH  (3)  =  5. 

an  array  of  length  NCH  which  specifies 
the  demagnification  factor  for  each  re¬ 
quested  channel. 

first  point  in  the  actual  sum,  predicted 
sum,  and  error  used  in  calculating  the 
power-spectra. 

an  integer  such  that  the  number  of  points 
for  computing  power-spectra  =  2n.  If 
N  -  0,  no  power  spectra  are  computed. 

an  integer  specifying  the  number  of 
times  to  smooth  the  power-spectra. 


Space  required:  24598  locations  used 


5.  Temporary  Storage  requirements ;  none 


6 .  Alarms :  none 


7>  Ihe  pr0gram  Prints  the  Program  name,  card 

mance  factor  S*  the  power-sPectr'a>  and  the  perfor- 

8‘  -E-rror  Printouts  and  stops:  if  any  of  the  following 
errors  occur,  the  type  Sf  error  is  printed  on 
the  program  is  stopped.  * 

a.  Seismogram  not  on  the  tape 

b.  100  many  channels  from  the  seismogram 

c.  Too  many  points  requested  from  the  seismogram 
•  Illegal  channel  requested  from  the  seismogram 

^ *  Input  and  output  tape  mountings: 

a)  Standard  run 

subset  tape  on  logical  unit  3.  Scratch  tapes 

“?lts  56>57, 2,  and  4.  Logical  unit  4 
is  the  plot  tape  and  logical  unit  2  is  the  save 
tape.  The  save  tape  contains,  in  binary  logical 
records,  the  following  information: 

Record  1:  ISEIS,  NCH,  NPTS,  ISPT,  IDETREND,  JUMP,  ISPAN 
LENGTH,  ICH(l)  to  ICH(NCH),  SF(1)  to  SF(NCH)  * 

EOF 

Record  2:  No • H?° jj t s ; in  the  correlation  matrix  (NCH2* 

Record  3:  Correlation  matrix 
EOF 

Record  4:  No.  of  points  in  the  filter  (NCH* LENGTH) 

Record  5:  Filter  matrix 
EOF 

Record  6:  NPTS  -  JUMP* (LENGTH- 1+ ISPAN) 

Record  7:  Sum  of  data  channels  (actual  array  sum) 

EOF 

Record  8:  NPTS  -  JUMP* ( LENGTH- 1+ I SPAN ) 

Record  9:  Predicted  array  sum 
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EOF 


Record  10:  NPTS  -  JUMP* (LENGTH- 1+ I SPAN) 

Record  11:  Error  (actual  array  sum  minus  predicted  array 
sum 

b)  Option  1 

Input  subset  tape  on  logical  unit  3.  Input  tape 
on  logical  unit  1,  (must  be  a  save  tape  in  the 
format  as  discussed  above).  Logical  units  2  and 
4,  56  and  57  are  scratch  tapes.  Unit  2  contains 
the  resulting  save  tape  in  format  described  above. 
Unit  4  is  the  plot  tape. 

c)  Option  2 

Input  subset  tape  on  logical  unit  3.  Input  tape 
on  logical  unit  1  (must  be  save  tape  in  the  format, 
as  discussed  above).  Scratch  on  tape  unit  4,  56, 
and  57.  Unit  4  is  the  plot  tape. 

10.  Input  and  output  card  formats:  The  format  for  card  1 
is  FORMAT  (915).  The  parameters  are  right  adjusted 
integers  and  are  in  the  following  fields: 

Columns  Parameters 


1-5 

ISEIS 

6-10 

NCH 

11-15 

NPTS 

16-20 

ISPT 

21-25 

I DETREND 

26-30 

JUMP 

31-35 

I  SPAN 

36-40 

LENGTH 

41-45 

I0PTI0N 

Card(s)  2  has  format  (1615).  The  parameters  are  right 
adjusted  integers  and  are  in  the  following  fields: 

Columns  Parameters 

ICH(l) 

ICH( 2 ) 

• 

•  • 

76-80  ICH (16 ) 

If  NCH  >16  a  second  card  must  be  used  with  ICH(17)  in  col¬ 
umns  (1-5),  ICH(18)  in  columns  (6-10),  etc. 


1-5 

6-10 
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Card(s)  3  has  format  (8F10.5).  The  parameters  are  float¬ 
ing  point,  and  the  numbers  should  be  punched  in  their  res¬ 
pective  fields  with  decimal  point. 

Columns  Parameters 


1-10 

SF(1 ) 

11-20 

SF(  2 ) 

• 

• 

• 

71-80 

• 

SF(  8 ) 

Repeated  cards  are  used  in  the  same  manner  as  Card  2 
if  NCH  >  8. 

Card  4  has  format  (315).  The  parameters  are  right  ad¬ 
justed  integers  and  are  in  the  following  fields: 


Columns 

I- 5 
6-10 

II- 15 


Parameters 

ISTART 

N 

NSM 


The  order  of  card  input  is  as  follows 
Standard  run:  Card  1 

Card  2 
Card  3 


parameters  describe  data  from 
which  filter  is  computed 


Option  1 


Card  1  "1  parameters  describe  data  to 
be  filtered 

Card(s )  2 
Card(s )  3 
Card  4 

Card  1  ]  parameters  describe  data  to 
be  filtered 

Card(s ) 2 
Card(s)3 
Card  4 
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Option  2: 


Card  1 


parameters  describe  data  to 
be  filtered 


Card(s )  2 
Card  ( s  )  3_ 

Card  4 

11.  Selective  jump  and  stop  settings:  none 

12.  Timing :  For  10  channels  of  data  with  2000  points  in 
each  channel. 

Standard  run:  30  minutes 

Option  1:  10  minutes 

Option  2:  5  minutes 

13.  Accuracy:  Single  precision 

1 4 .  Caution  to  users : 

NCH  <  16 

NCHCof  data  for  filter  computation)=NCH(of  data  being 
filtered) 

LENGTH  <40 
ISPAN  £  10 
NPTS  <  4000 
N  <  10 

15.  Equipment  configuration: 

Standard  run:7/9  COOP ,XXXXX, I/BY  1/3/S/  2/4/ 56/ 57/E/ 5= 56 , 60 , 5000 

Option  1:  7/9  COOP ,XXXXX,I/  1/3/S/  2/ 4/ 56/ 57/E/ 5= 56 , 60 , 5000 

Option  2:  7/9/  COOP ,XXXXX,I/  1/ 3/ S/BY2/4/ 56/ 57/E/ 5= 56 , 60 , 5000 

16 .  References : 

(1)  Teledyne  Internal  Memorandum  from  E.  A.  Flinn 
and  J.  Claerbout  to  R.  Van  Nostrand,  dated  16 
September  1965.  Subject  "Some  topics  in  digital 
filtering  theory  and  application". 

(2)  McCowan,  D.W.,  Teledyne  Technical  Memorandum  No. 
8-66,  dated  16  December  1966,  Title:  "Finite 
Fourier  Transform  Theory  and  Its  Application  to 
the  Computation  of  Convolutions,  Correlations, 
and  Spectra". 


-57- 


D.  METHOD 


The  steps  below  show  the  general  method  used  for  com¬ 
putation  . 

(1) .  Compute  correlation  matrix  of  correlation  functions: 

N  N  M-l 

R(i,j,L)  =  l  l  l  X.(K)  X.CK+L) 

i=l  j=l  k=0 

where  M  is  the  number  of  points  per  channel,  N  is  the 
number  of  channels,  and  L  is  the  lag.  See  reference  (1). 

(2) .  Compute  prediction  filter  using  the  Levinson-Robinson 
recursvion  method. 

C  3 ) .  Filter  data 

(4) .  Plot  direct  sum  of  data  channels  (actual  sum),  and 
the  predicted  sum.  Also  plot  error,  where  error  =  actual 
sum  -  predicted  sum. 

(5)  Compute  power-spectra  of  the  three  traces  above.  See 
reference  (2). 

(6)  Plot  power-spectra  traces  and  a  performance  trace 
where  performance  =  1  -  error  spectra/actual  spectra. 
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SEISMIC  DATA  LABORATORY 
Alexandria , Virginia 
Digital  Computing  Section 


A.  IDENTIFICATION 

Title:  Spectral  Matrix  Averaging  Program 

COOP  Identification:  G635  SPECAVE 
Category:  Time  Series  Analysis 

Programmer:  D.  W.  McCowan 

Date:  October  iSo7 

B.  PURPOSE 

This  program  computes  and  averages  spectral  matrices  for 
multichannel  seismic  array  data*  The  output  of  this  program 
is  a  tape  in  the  same  format  as  the  program  HEFALUMP  save  tape. 
Using  that  tape  as  input  to  the  HEFALUMP  program  allows  the 
user  to  compute  multichannel  filters  from  an  ensemble  of  noise 
spectra.  This  program  is  also  designed  to  add  to  or  update  an 
ensemble  spectral  matrix  on  tape  from  a  previous  run. 

C.  USAGE 


!•.  Operational  procedure:  This  is  a  Fortran-63  main  pro¬ 
gram  with  the  following  subroutines:  MAKESUB,  DETRND,  COOLER, 
SMOOTH,  and  XMNG .  In  addition,  the  following  utility  sub¬ 
routines  are  assumed  to  be  on  the  system  tape:  ERASE,  COOL, 
DISCG3 ,  ana  SKIPREC. 

2.  Parameters:  The  user  is  referred  to  the  write-up  of 
program  HEFALUMP  for  a  clearer  explanation  of  some  of  these 
parameters . 

a)  NCASES  -  The  number  of  spectral  matrices  to  com¬ 

pute  and  add  to  the  ensemble  in  this  run. 

b)  ISW  -  The  run  designation  switch. 

ISW=0  only  the  NCASES  spectral  matrices 
computed  this  time  will  be  on  the  save 
tape . 

ISW^O  the  NCASES  spectral  matrices  com¬ 
puted  this  time  will  be  added  to  the 
ensemble  input  on  a  save  tape  from  a 
previous  run  mounted  on  unit  4  and  writ¬ 
ten  on  the  save  tape. 
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c)  ISEIS 

d)  NCH 

e)  NPTS 

f)  ISPT 

g)  IDT 

h)  L 

i)  IPAUSE 


j)  ICH(I) 

k)  X(I) 

l)  Y (I ) 

m)  srci) 


The  seismogram  number  of  this  sample  of 
data, 


The  number  of  channels  in  this  sample 
of  data.  This  must  be  the  same  for  each 
sample.  NCH  <_  13. 

The  number  of  points  in  this  sample  of 
data.  NPTS  must  be  a  power  of  two. 

NPTS  <_  409  6. 

The  first  point. 

The  detrend  switch  for  this  sample. 

IDT  =  -1  for  no  detrending. 

IDT  =  0  to  remove  the  mean. 

IDT  =  1  to  remove  the  mean  and  linear 
trend. 

The  number  of  points  desired  in  the 
spectral  estimates.  L  must  be  a  power 
of  two  plus  1  and  must  be  the  same  for 
each  sample.  L  £  129. 

A  pause  switch. 

IPAUSE  =  0  has  no  effect. 

IPAUSE  =  0  pauses  with  the  index  number 
of  the  sample  being  computed  in  the  A 
register  to  allow  a  new  input  subset 
tape  to  be  mounted. 

The  channel  index  of  the  Ith  data  channel 
of  seismogram  ISEIS.  ICH(I)  =  4  indicates 
the  Ith  data  channel  is  the  fourth  on  the 
input  tape,  etc. 

The  X_ coordinate  of  the  Ith  data  channel. 
This  is  optional  input. 

The  Y . coordinate  of  the  Ith  data  channel. 
This  is  optional  input. 

The  demagnification  factor  for  the  Ith 
data  channel.  This  number  is  divided 
into  every  point  of  the  Ith  data  channel. 
If  it  is  0.0  it  is  set  to  1.0. 


3.  Space  required:  32K  locations 
4-  Temporary  storage;  DISC 

5.  Alarms :  Various  alarms  for  illegal  data  are  printed 
out . 
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6. 


Error  returns: 


none 


7*  Error  stops:  A  stop  occurs  after  each  alarm. 

8*  Taye  Mountings:  Tape  unit  1  is  the  input  subset  tape 
which  can  be  changed  during  execution  by  using  the 
pause  feature .  Tape  unit  2  is  the  output  save  tape 
containing  the  sum  of  all  computed  spectra  and  any 
input  from  a  previous  run.  Tape  3  is  a  scratch  tape. 
Tape  4  is  an  input  save  tape  from  a  previous  run  of 
this  program,  if  so  desired.  If  not,  it  should  be 
neglected  with  the  BY  feature  of  the  monitor. 

9.  Input  and  output  formats:  Input  cards: 

Parameters  (a)  and  (b)  appear  on  the  first  card  as 
FORMAT  (2110).  There  follows  NCASES  data  decks  of 
the  form:  card  1  containing  parameters  (c)  through 
(i)  as  FORMAT  (110,615),  and  the  next  cards  contain¬ 
ing  parameters  (j)  through  (m)  as  FORMAT  (3(13,  2F8.2, 
F6  .  2 )  ) . 


Save  tape:  The  save  tape  consists  of  the  two  logical 
records  of  binary  data.  The  first  record  is  a  127 
word  binary  label  written  as  follows: 


LAB ( 1 )  =  ISEIS 
LAB( 2 )  =  NCH 
LAB ( 3 )  =  NPTS 
LAB (99)  =  NCASES 


LAB ( 100 )  =  IDT 
LAB (101)  =  L 
LAB(58  — 58+NCH-l)  =  X 
LAB(79  —  79+NCH-l)  =  Y 


The  values  of  NCH,  NPTS,  IDT,  L,  X  and  Y  that  appear 
in  the  label  are  those  of  the  last  data  sample.  NCASES 
is  always  updated  in  this  program  and  thus  is  an  ac¬ 
curate  number  of  the  separate  spectral  matrices  in  the 
ensemble.  ISEIS,  in  true  fashion,  is  the  average  of 
the  values  read  in  for  the  samples. 


The  second  record  is  a  21801  word  array  with  the 
averaged  ensemble  spectral  matrix  stored  close  packed 
and  written  as  follows: 


sij(")  * 


S(k,j,w),  ( (S(i, j ,w) ,  i  =  1,  NCH),  j  =1 ,  NCH),  w=l,L) 


10. 

Selective 

jumps : 

none 

11. 

Timing: 

See  write 

s-up  of  subroutine  XMNG 

12. 

Accuracy : 

Single 

precision 

J 
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3.  Cautions  to  user:  none 

4-  Equipment  configuration:  Standard  COOP  with  RELOCOM. 
5 .  References : 

Write-ups  of  program  HEFALUMP  and  subroutine  XMNG 
by  D.  W.  McCowan  and  J.  Jih  respectively. 

METHOD 


See  references. 


SEISMIC  DATA  LABORATORY 


Alexandria , Virginia 
Digital  Computer  Section 


A.  IDENTIFICATION 

Title:  Spatial  interpolation 

COOP  Identification:  G625  TWX 
Category:  Time  Series  Analysis 

Programmer:  E.  A.  Flinn 

Date:  13  February  1967 

B.  PURPOSE 

To  compute  a  least-mean-square-error  filter  which  inter¬ 
polates  one  channel  of  an  array  from  certain  other  channels 
in  the  array,  apply  the  filter  to  construct  an  interpolated 
estimate,  construct  the  error  time  function,  and  display  the 
power  spectrum  of  the  error. 

C.  USAGE 


1.  Operational  procedure:  This  is  a  main  program  written 
in  Fortran-63.  Relocom  and  Executer  are  both  required. 

2.  Parameters :  The  following  control  parameters  are  in¬ 
put  from  cards  (see  C.9  below): 

a)  NSEIS  -  Seismogram  number  containing  data  to  be 

filtered.  The  data  is  read  in  from  a 
subset  tape,  but  the  specified  seismo¬ 
gram  does  not  have  to  be  the  first 
seismogram  on  the  tape.  An  error  stop 
occurs  if  the  seismogram  is  not  on  the 
tape  at  all. 


b) 

NCH 

' 

Number  of  channels  to  be  re-subset  from 
the  input  subset  tape. 

c) 

11 

- 

Start  point  for  re-subset. 

d) 

12 

- 

Number  of  points  to  re-subset, 
stop  occurs  if  11+12  is  greater 

Error 

than 

the  third  word  of  the  subset  seismogram 
label . 
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e)  CHN 


f )  NSM 


g )  LAGS 


h)  IPLOTR 


i)  I CHN 


j)  ICOR 


Hollerith  designator  for  the  channel  to 
be  interpolated.  Error  stop  occurs  if 
label  of  input  subset  tape  does  not  con¬ 
tain  this  designator. 

Number  of  times  to  decimate  the  filtered 
data  by  two,  in  constructing  the  spectra. 
Only  the  first  N  points  of  the  filtered 
trace  are  used  in  forming  the  spectra, 
where  N  is  the  largest  power  of  2  less 
than  the  number  of  data  points.  For  ex¬ 
ample,  1035  input  points  are  reduced  to 
1016  by  filtering  with  a  20-point  filter. 
The  program  takes  the  first  512  data 
points  in  forming  the  spectra.  Suppose 
NSM  =  4:  then  the  512  points  are  reduced 
to  512/4  =  32  spectral  points  which  are 
plotted. 

Number  of  elements  in  the  filter  and 
number  of  lags  in  the  correlation  matrix. 
Maximum  is  50. 

Plot  control  switch  for  the  correlation 
matrix;  input  zero,  correlation  matrix 
is  plotted;  input  non-zero,  omit  plot. 

Channel  number  to  interpolate,  counted 
after  re-subsetting.  See  definition  of 
ISUB  below.  For  example,  suppose  there 
are  5  channels  in  the  input  tape;  these 
are  taken  to  be  numbered  1,  2,  3,  4,  5. 

If  we  wish  to  interpolate  channel  3  from 
channels  1  and  we  make  NCH  =  3  (see 
above),  ISUB(1,  2,  3)  =  (1,  3,  4)  (.see 
below),  and  put  ICHN  =  2,  since  the 
channel  to  be  interpolated  is  the  second 
channel  on  the  re-subset  tape. 

Switch  for  computing  correlation  matrix 
and  filter  or  reading  them  in  off  tape. 
ICOR  =  0  causes  correlation  matrix  and 
filter  to  be  computed  and  saved  on  a 
binary  output  tape  on  unit  9.  (see  C.9 
below) . 

ICOR  =  1  makes  the  program  read  in  the 
correlation  matrix  and  filter  coefficients 
from  binary  tape.  An  error  stop  occurs 
if  NCH  is  not  equal  to  the  second  word  of 
the  label  of  the  tape  on  9,  and  if  LAGS 
is  not  equal  to  the  third  word  of  the 
label  of  the  tape  on  9. 
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The  filter  is  recomputed  if  CHN  is  not 
equal  to  the  37th  word  of  the  label  of 
the  tape  on  9. 


In  any  case,  the  tape  on  9  is  re-wound 
and  rewritten  (see  C.9  below) 

This  allows  the  plots  to  be  less  than 
full-length, ie,  I2-I1-LAGS  points.  Non¬ 
zero,  only  IPLOTX  points  are  plotted. 

An  error  stop  occurs  if  IPLOTX  exceeds 
I  2- II- LAGS. 

An  integer  array  of  channel  numbers  to 
re-subset.  These  must  be  arranged  in 
increasing  numerical  order.  If  all  chan¬ 
nels  are  to  be  re-subset,  ISUB  may  be 
left  blank. 

An  array  of  channel  magnifications,  ar¬ 
ranged  in  the  same  order  as  ISUB.  Each 
re-subset  channel  is  divided  by  the 
corresponding  magnification.  Any  element 
of  SF  left  blank  is  set  to  1.0  by  the 
program. 

3.  Space  required:  32K  using  Relocom  and  Executer. 

4.  Temporary  locations:  none 

5.  Alarms  and  special  printouts :  none 

6.  Error  returns;  not  applicable. 

7.  Error  stops:  An  error  message  is  printed  and  the  pro¬ 
gram  stops  under  the  following  circumstances. 

a.  Requested  seismogram  not  on  input  subset  tape. 

b.  Number  of  channels  requested  greater  than  number 

of  channels  in  requested  seismogram  or  greater  than 
16. 

c.  Number  of  points  requested  greater  than  number  of 
points  in  requested  seismogram. 

d.  Label  of  requested  seismogram  does  not  contain 
channel  with  designator  CHN  (see  C.2  above). 


k)  IPLOTX 


1)  ISUB 


m)  SF 


e.  Number  of  lags  input  exceeds  50. 


f .  When  correlation  matrix  is  read  in  from  input 
tape,  NCH  different  from  second  word  of  label 
or  LAGS  different  from  third  word  of  label. 

g.  ISUB  not  in  increasing  numerical  order. 


8 .  Tape  mountings: 


Input  tapes 


1  Input  subset  tape. 

7  Card  input. 

9  Correlation  matrix  and  filter 
(optional ) . 


Output  tapes  - 

2  Plot. 

6  Printer. 

9  Correlation  matrix  and  fi.  ter. 


Scratch  tapes:  3,4,5 


The  COOP  card  arrangement  for  using  on-line  card  reader  and 
printer  is:  I/l/0/2/S/3/4/5/9/E/7=50/6=51/8=54. 


9. 

Input  and 

Card  1: 

output  formats:  Input  cards  - 

Format  (415,  2X,  A3,  615) 

Columns 

Data 

Explanation 

1-5 

NSEIS 

Seismogram  desired. 

6-10 

NCH 

Number  of  channels  to  be  re-subset, 
(maximum  is  16). 

11-15 

11 

Starting  point  for  re-subset. 

16-20 

12 

Number  of  points  to  be  re-subset. 

23-25 

CHN 

Hollerith  designator  for  channel  to 
be  interpolated. 

26-30 

NSM 

Number  of  times  to  decimate  spectra 
by  2  before  plotting. 

31-  3  5 

LAGS 

Number  of  elements  in  filter  (maximum 
is  50 ) . 

36-40 

IPLOTR 

=  0  plot  correlation  matrix; 
i  0  omit  plot. 

41-45 

I  CHN 

Order  number  of  CHN  on  re-subset  tape 

46-50 

ICOR 

=0  compute  correlation  matrix  and 
filter,  and  write  on  tape  9. 
i  0,  read  in  correlation  matrix  and 
filter  from  tape  9  (filter  is  recom¬ 
puted  if  CHN  differs  from  word  37  of 
label  of  tape  9 ) . 

51-55 

IPLOTX 

=  0  plot  channel  CHN,  interpolated 
estimate  of  CHN, and  error,  all  full- 
length; 
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negative,  omit  above  plots, 
positive,  plot  only  IPLOTX  points. 

Card  2:  Format  (2013) 

This  card  contains  ISUB,  an  integer  array  of  channel  numbers 
to  be  re-subset.  An  error  stop  occurs  if  these  are  not  in  in¬ 
creasing  numerical  order.  For  example,  to  re-subset  channels 
1,  3,  and  4,  punch  1  in  column  3,  3  in  column  6,  and  4  in  col¬ 
umn  9 . 

Card  3  and  4:  Format  (12F5.Q) 

These  cards  contain  SF,  the  array  of  channel  magnificat¬ 
ions  for  the  re-subset  tape.  The  magnifications  must  corres¬ 
pond  to  the  elements  of  ISUB.  In  the  above  example,  columns 
1-15  of  card  3  would  contain  the  magnifications  for  channels 
1,  3,  and  4;  card  4  would  be  blank.  Any  magnifications  left 
blank  will  be  set  to  1.0  by  the  program.  Thus  if  no  demagni¬ 
fication  is  required,  cards  3  and  4  may  both  be  left  blank. 

Input  Tapes  - 

1:  A  data  tape  in  subset  format,  containing  seismogram 

number  NSEIS . 

9:  See  output  formats  below. 

Scratch  tape- 

5:  A  re-subset  version  of  the  input  tape  1,  containing 

only  NCH  channels  as  specified  by  ISUB,  i.e.,  the  label  is 
the  same  as  the  label  on  tape  1,  with  the  following  except¬ 
ions:  word  2  =  NCH;  word  3  =  1 2-1 1+1;  the  first  NCH  words 

of  the  channel  designator  group  (beginning  at  word  8)  are 
the  designators  of  the  NCH  channels  as  specified  in  ISUB. 

This  tape  is  later  overwritten  by  the  program. 

Output  Tapes  - 

2.  Plot  tape: 

a.  Correlation  matrix  (if  IPLOTR  =  0),  displayed 

in  standard  matrix  form:  half  the  channel  auto¬ 
correlations  down  the  diagonal,  the  right  half 
of  the  crosscorrelations  above  the  diagonal,  and 
the  (reversed)  left  halves  of  the  cross-correlations 
below  the  diagonal. 

b.  Input  channel  CHN. 

c.  Interpolated  estimate  of  CHN. 

d.  Error  -  (b)  minus  (c). 

(b,c,d  are  all  to  the  same  scale  factor,  with  range  h  inch. 

The  maximum  value  in  CHN  is  written  in  BCD  on  the  plot) 
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fh  .  ,  e*  E£ror’’  with  its  own  scale  factor. 

unless "iPLOTX^are "non- zero' ' £ ^ • 

IPLOTX  points)  9  whlch  case  the  plots  contain 

f.  Spectrum  of  CHN;  maximum  is  identified 
is  identified^  Um  °f  interP°luted  estimate  of  ON;  maximum 

identified.'  SpeCtrum  of  error  (*»«  d  above);  maximum  is 

plotted  linearly?  ^he^eft^nd^f^he^jl”?1*  range  1!l  inches . 
zero  frequency  the  SLC; I  of*he  Plot  corresponds  to 

the  sampling  given  in  the  label°  ^Th  folddng  frecluency  -  half 
plots,  f.  g  h.  i  i  K  =  m  5  v  "UI!lber  of  Points  in 
of  2  less  than  I2-I1-LAGS.)  *  where  N  1S  the  largest  power 
i.  Performance  factor,  defined  =,=,  i  n 
of  interpolated  spectrum  to  input  spectrum, 


9: 


Correlation  matrix  and  interpolation  filter: 
Record  1 


Record  2; 


Record  3: 


Printer  - 


Label,  identical  to  label  on  scratch  tape  5 
above,  except  that  word  37  =  CHN. 

Correlation  matrix,  written  with  the  state 
ment  WRITE  (9)  «RU,J,K),  1=1,  NCH)  ,J=l?NCH) , 

Interpolation  filter,  written  wi+h  thp  + 
ment  WRITE  (9)  CCSCl.J),  1=1, NCH), 


channeircHN?  ^hfinfeipo^tfre^tima^r^f^SN8  ^d^he'"^1 

23  STSSf^S: 

10 ■  Jump  settings:  none 

11-  TimeRequired:  depends  on  the  number  of  channels 
and  the  length  of  the  filter.  Time  is  dominoed 
„ppf  ,bJ  computation  of  the  correlation  matrix  and 
second  by  computation  of  the  filter  =  i  d 

gjpafiS^f— 

filter,  and  3  minutes  for  everything  else  s  IS* 
nels  and  20  lags  requires  15  minufef.  °han‘ 

12.  Accuracy:  single  precision 


13 


14, 


Cautions  to  user:  Refer  to  error  messages  above 
Equipment  configuration;  standard  COOP 
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15.  References : 


SDL  Report  No.  168  "Finite  Fourier  Transform  Theory 
and  its  Application  to  the  Computation  of  Convolutions, 
Correlations,  and  Spectra",  by  D.  W.  McCowan,  15 
December  1966. 
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SEISMIC  DATA  LABORATORY 


Alexandria,  Virginia 
Digital  Computer  Section 

A.  IDENTIFICATION 
Title:  VFKSPTRUM 

COOP  Identification:  G629  VFKSPTRUM 
Category:  Time  Series  Analysis 
Programmer:  J.  Jih 
Date:  2  June  1967 

B.  PURPOSE 


To  compute  and  display  the  frequency-wave 
spectra  of  seismic  noise  along  with  a  response 
the  corresponding  vertical  array. 


function 


C.  USAGE 


lm  Operational  procedure:  This  is  a  Fortran- 6 3  main 
program  with  the  following  subroutines:  FKMATRIX  SMOOTH 

™D,and  CNTUR6.  In  a"S,S!Se™;,- 
lowing  utility  subroutines  are  assumed  to  be  on  the  svstem 
tape:  SKIPREC,  ERASE,  DISC63,  and  COOL.  7 


2*  Parameters :  This  program  processes  job 
sequence,  reading  m  the  required  data  cards  each 
3 oh  requires  the  following  data  cards: 


requests  in 
time.  Each 


A. 


FORMAT**^  6I10^2F8?5^  ^^followinB  Parameters  as 


a. 

b. 


ISM 

N 


-  the  seismogram  number  of  the  desired  data 

-  the  number  of  data  channels  to  be  used. 
1<N<  24 . 


c .  LX 


d.  ISPT 


the  number  of  digital  points  to  be  used 
from  each  data  channel.  LX  must  be  a 
power  of  two.  (32  <  LX  <  4096).  If  LX 
>  4096,  it  will  be  truncated  to  4096,  If 
LX  <  32,  it  will  be  set  to  32. 

the  first  point  of  the  desired  data. 

(1  <  ISPT). 
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e.  IDT  -  the  detrend  switch. 

=  0  or  blank  to  remove  only  the  mean 
=  1  to  remove  the  mean  and  linear  trend 
=  -1  for  no  detrending. 

f.  NK  -  the  number  of  K,  vertical  wave  number, 

values  to  compute.  (1  <_  NK  <_  99). 


g.  SK  -  the  lowest  K  value  to  compute 

h.  DK  -  the  increment  between  successive  K 

values . 


i.  IDEC  -  a  decimation  parameter.  AN  F-K  contour 

between  DC  and  the  (folding  frequency) 
*2**(1-IDEC)  will  be  computed  and  dis¬ 
played.  If  IDEC  <  1,  it  will  be  set 
to  1.  2** (IDEC+4 )  <  LX. 

j  .  IRESP  -  a  response  switch. 

=  0  or  blank  for  no  array,  response 
i  0  for  a  plot  of  the  array  response. 

B.  The  next  cards  list  the  following  parameters  as  FORMAT 
(4(12, F9. 3,  F9 . 3  ) ) : 


a.  ICH(I) 

b.  Y(I) 


the  position  of  the  Ith  desired  data 
channel  on  the  input  seismogram.  (1  <_  I  <_  N). 

the  depth  of  the  Ith  desired  data  chan¬ 
nel.  (1  <  I  <  N): 


c.  F(I)  -  scale  factor  for  the  Ith  channel,  which 

is  divided  into  the  data,  (1  <  I  <  N). 
It  will  be  set  to  1,  if  it  is  zero  or 
blank . 

3 .  Space  required :  20480 


4.  Temporary  space :  DISC 

5.  Alarms :  See  write-up  of  "MAKESUB" . 


6.  Error  returns :  none 

7.  Error  stops:  There  are  two  error  stops: 

a.  The  control  data  cards  are  not  enough. 

b.  N,  NK  or  DK  is  illegally  requested. 
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8.  Tape  mountings:  An  input  tape  of  subset  seismograms 


must  be  on  logical  Unit  1.  A  scratch 
ed  and  must  be  on  logical  Unit  2. 

tape  is  need- 

9  . 

Selective  jumps:  none 

10. 

Timing:  Example  N  =  5,  LX  =  4096, 

IRESP  =  1,  case  takes  7  minutes  of 

NK 

CDC 

=  21,  IDEC  =  4, 
1604B  time. 

11. 

Accuracy:  Single  precision 

12. 

Caution  to  user:  none 

13.  Equipment  configuration:  Standard  COOP 

14.  References :  See  write-up  of  "FKSPTRUM" . 


D.  METHOD 


The  frequency-wave  number  power  spectra  are  computed  from 
the  spectral  matrix  elements  by  the  following  relation: 


P(f,k)  =  l 
i 


l  Sij(f)t 

j 
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